{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fffbbaec",
   "metadata": {},
   "source": [
    "# 利用 Product Formula 模拟时间演化\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6f9a147",
   "metadata": {},
   "source": [
    "## 概述\n",
    "\n",
    "量子力学中系统的能量由哈密顿量算符 $H$ 描述，求解给定系统的哈密顿量的全部或者部分性质，构成了凝聚态物理、计算化学和高能物理等一系列学科的核心问题。然而，由于系统的自由度随系统增大呈指数级增加，导致一般情况下无法利用经典计算机有效模拟量子系统——即便用光全世界的内存和硬盘也不能直接储存一个仅几百量子比特的系统的态矢量。与经典计算机不同，量子计算机由于其所有的操作都是直接作用在同样为指数级别的态空间上，因而在模拟一个量子系统上量子计算机具有经典计算机无法比拟的优势。实际上，设计一个可控的量子系统来高效模拟自然界中的量子系统，正是费曼在上世纪 80 年代提出量子计算这一概念时的最初想法：\n",
    " \n",
    "> _\"Nature isn't classical, dammit, and if you want to make a simulation of nature, you'd better make it quantum mechanical, and by golly it's a wonderful problem, because it doesn't look so easy.\"_\n",
    ">\n",
    "> --- \"Simulating physics with computers\", 1982, Richard P. Feynman [1]\n",
    "\n",
    "通用量子计算机以及一系列量子模拟器的发展令费曼的设想有了实现的可能。在通用量子计算机上进行数字量子模拟（digital quantum simulation）—— 利用量子门构造量子线路从而实现量子模拟，由于该方法具有较高的可拓展性和通用性，因而被认为是最有潜力的技术路线。\n",
    "\n",
    "本教程讲述了如何利用 Paddle Quantum 模拟量子系统的时间演化过程，主要分为以下三个部分：\n",
    "1. 如何在 Paddle Quantum 中创建并操作一个 `Hamiltonian` 对象\n",
    "2. 如何利用 `construct_trotter_circuit` 来构建时间演化电路\n",
    "3. Suzuki product formula 方法的原理，以及如何进一步搭建任意阶的 Trotter-Suzuki 电路\n",
    "\n",
    "\n",
    "## 定义系统哈密顿量\n",
    "\n",
    "在介绍如何搭建时间演化电路之前，读者可以先熟悉一下 Paddle Quantum 中的哈密顿量。目前用户需要通过定义一个列表的形式来创建哈密顿量，列表中的元素应为哈密顿量中的每一项的系数与其泡利算符。作为教程，我们首先考虑一个简单的哈密顿量：\n",
    "\n",
    "$$\n",
    "H = Z \\otimes Z\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "该哈密顿量描述了两个量子比特之间的一种简单相互作用：当两个量子比特同时处于 $|0\\rangle$ 态或者 $|1\\rangle$ 态时，系统的能量为 $+1$；相反，当两个量子比特的态不同时，系统的能量为 $-1$。\n",
    "\n",
    "接下来，用户可以在 Paddle Quantum 中创建这一哈密顿量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1a17d7d8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 Z0, Z1\n"
     ]
    }
   ],
   "source": [
    "from paddle_quantum.utils import Hamiltonian\n",
    "\n",
    "h = Hamiltonian([[1, 'Z0, Z1']])\n",
    "print(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64ef6a63",
   "metadata": {},
   "source": [
    "目前，Paddle Quantum 中的哈密顿量类 `Hamiltonian` 支持自动合并同类项，加减法、索引以及拆分等操作："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ff08a2a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 Z0, Z1\n"
     ]
    }
   ],
   "source": [
    "h = Hamiltonian([[0.5, 'Z0, Z1'], [0.5, 'Z1, Z0']], compress=True)\n",
    "print(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fec891a5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h + h: 2.0 Z0, Z1\n",
      "h * 2: 2.0 Z0, Z1\n",
      "h: 1.0 Z0, Z1\n"
     ]
    }
   ],
   "source": [
    "print('h + h:', h + h)\n",
    "print('h * 2:', h * 2)\n",
    "print('h:', h[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2be41b1",
   "metadata": {},
   "source": [
    "同时，内置的 `decompose_pauli_words()` 和 `decompose_with_sites()` 方法可以将哈密顿量分解为更加方便处理的形式："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "40fcc0b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pauli words 分解： ([1.0], ['ZZ'])\n",
      "Pauli with sites 分解： ([1.0], ['ZZ'], [[0, 1]])\n"
     ]
    }
   ],
   "source": [
    "print('Pauli words 分解：', h.decompose_pauli_words())\n",
    "print('Pauli with sites 分解：', h.decompose_with_sites())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95945ddd",
   "metadata": {},
   "source": [
    "除此之外，`construct_h_matrix()` 还可以创建其在泡利 $Z$ 基底下的矩阵形式："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "497ff8fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.+0.j,  0.+0.j,  0.+0.j,  0.+0.j],\n",
       "        [ 0.+0.j, -1.+0.j,  0.+0.j,  0.+0.j],\n",
       "        [ 0.+0.j,  0.+0.j, -1.+0.j,  0.+0.j],\n",
       "        [ 0.+0.j,  0.+0.j,  0.+0.j,  1.+0.j]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h.construct_h_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "719696dd",
   "metadata": {},
   "source": [
    "## 模拟时间演化\n",
    "\n",
    "根据量子力学的基本公理，在确定了一个系统的哈密顿量之后，该系统随时间演化的过程可以由如下方程描述\n",
    "\n",
    "$$\n",
    "i \\hbar \\frac{\\partial}{\\partial t} | \\psi \\rangle = H | \\psi \\rangle,\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "$\\hbar$ 为约化普朗克常数。该方程正是著名的薛定谔方程。因此，对于一个不含时的哈密顿量，系统的时间演化方程可以写为\n",
    "\n",
    "$$\n",
    "|\\psi(t) \\rangle = U(t) | \\psi (0) \\rangle, ~ U(t) = e^{- i H t}.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "这里我们取自然单位 $\\hbar=1$，$U(t)$ 为时间演化算符。利用量子线路来模拟时间演化过程的核心思想是利用量子电路构建出的酉变换模拟和近似该时间演化算符。在量桨中，我们提供了 `construct_trotter_circuit(circuit, Hamiltonian)` 函数来构建对应某一个哈密顿量的时间演化电路。下面，就让我们用刚刚创建的哈密顿量来测试一下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3449029b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--*-----------------*--\n",
      "  |                 |  \n",
      "--x----Rz(2.000)----x--\n",
      "                       \n"
     ]
    }
   ],
   "source": [
    "from paddle_quantum.trotter import construct_trotter_circuit\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "\n",
    "cir = UAnsatz(2)\n",
    "construct_trotter_circuit(cir, h, tau=1, steps=1) \n",
    "print(cir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3dc6c7f",
   "metadata": {},
   "source": [
    "可以看到，我们已经创建了针对 `h` 的基本量子模拟电路，它可以根据输入的 `tau` 来模拟任意时间长度的时间演化。\n",
    "\n",
    "此时，如果检查该电路对应的酉变换的矩阵形式，会发现它与时间演化算符 $e^{-iHt}$ 完全一致。这里，我们先使用 `gate_fidelity` 来计算量子电路的酉矩阵和时间演化算符的酉矩阵之间的保真度。保真度越接近 1 时代表两个酉矩阵代表的演化过程越相似。在下文，我们还会引入更加严格的误差定义。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8d312985",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "电路的酉算符和时间演化算符之间的保真度为： 1.000\n"
     ]
    }
   ],
   "source": [
    "from scipy import linalg\n",
    "from paddle_quantum.utils import gate_fidelity\n",
    "\n",
    "# 计算 e^{-iHt} 和电路的酉矩阵之间的保真度\n",
    "print('电路的酉算符和时间演化算符之间的保真度为： %.3f' \n",
    "      % gate_fidelity(cir.U.numpy(), linalg.expm(-1 * 1j * h.construct_h_matrix())))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa34ba8a",
   "metadata": {},
   "source": [
    "实际上，这是因为泡利算符组成的张量积对应的任意角度旋转算符都可以被很好的分解为门电路的形式。比如在这个例子中我们只需要改变电路中 Rz 门的角度，就可以实现任意 $e^{-i Z\\otimes Z t}$ 演化算符的模拟。那么，这是否意味着可以用类似的电路形式来对任意的写成泡利形式的哈密顿量进行完美的模拟呢？答案是否定的。考虑一个稍微复杂一些的哈密顿量：\n",
    "\n",
    "$$\n",
    "H = Z \\otimes Z + X \\otimes I + I \\otimes X.\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "同样地，让我们用 `construct_trotter_circuit` 来构建它所对应的时间演化电路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "40bcf1e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--*-----------------*----Rx(2.000)--\n",
      "  |                 |               \n",
      "--x----Rz(2.000)----x----Rx(2.000)--\n",
      "                                    \n",
      "电路的酉算符和时间演化算符之间的保真度为： 0.681\n"
     ]
    }
   ],
   "source": [
    "h_2 = Hamiltonian([[1, 'Z0, Z1'], [1, 'X0'], [1, 'X1']]) # 不需要写出单位算符\n",
    "cir = UAnsatz(2)\n",
    "construct_trotter_circuit(cir, h_2, tau=1, steps=1)\n",
    "print(cir)\n",
    "print('电路的酉算符和时间演化算符之间的保真度为： %.3f' \n",
    "      % gate_fidelity(cir.U.numpy(), linalg.expm(-1 * 1j * h_2.construct_h_matrix())))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f0a4ed",
   "metadata": {},
   "source": [
    "可以看到，此时电路对应的酉矩阵和时间演化矩阵之间的保真度 $<1$，说明该电路无法正确地模拟系统的时间演化过程。\n",
    "\n",
    "其实，对于我们此时构建的电路而言，它对应的酉变换为 $e^{-iZ\\otimes Z t} e^{-i (X\\otimes I + I\\otimes X)t}$，而时间演化算符为 $e^{-iZ\\otimes Z t - i(X\\otimes I + I\\otimes X)t}$。而对于一个量子系统而言，当两个算符不对易，即 $[A, B] \\neq 0$ 时，$e^{A+B} \\neq e^A e^B$。在这里，因为泡利算符 $X$ 和 $Z$ 之间并不对易，所以电路对应的酉变换并不等于正确的时间演化算符。\n",
    "\n",
    "除了利用保真度来描述量子电路和希望模拟的时间演化算符之间的相似性之外，我们还可以定义如下的误差 $\\epsilon$ 来刻画模拟演化电路和演化算符之间的距离，\n",
    "\n",
    "$$\n",
    "\\epsilon(U) = \\Vert e^{-iHt} - U \\Vert,\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "其中 $\\Vert \\cdot \\Vert$ 表示取最大的本征（奇异）值的模。这样的定义比起保真度而言更好的描述了量子态在不同的演化算符作用下产生的偏差，是一种更加严谨的定义模拟时间演化误差的方式。我们在下文中也会多次用到该误差。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a0f9969",
   "metadata": {},
   "source": [
    "### Product formula 与 Suzuki 分解\n",
    "\n",
    "Seth Lloyd 在 1996 年的文章中指出，可以将一整段的演化时间 $t$ 拆分为 $r$ 份较短的“时间块”来减小模拟时间演化的误差 [2]。考虑一个更一般的哈密顿量形式 $H = \\sum_{k=1}^{L} h_k$，其中 $h_k$ 是作用在一部分系统上的子哈密顿量。通过泰勒展开，不难发现其模拟误差是一个二阶项，即\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\prod_{k=1}^{L} e^{-i h_k t} + O(t^2).\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "那么，我们令 $\\tau = t/r$，并考虑演化算符 $\\left(e^{-iH \\tau}\\right)^r$，其演化误差为\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\left(e^{-iH \\tau}\\right)^r = \\left(\\prod_{k=1}^{L} e^{-i h_k \\tau} + O(\\tau^2) \\right)^r = \\left(\\prod_{k=1}^{L} e^{-i h_k \\tau} \\right)^r + O\\left(\\frac{t^2}{r}\\right).\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "上式告诉我们，只要将一整段演化时间拆为足够多的“片段”，就可以达到任意高的模拟精度，这就是 product formula 的基本思想。不过，(7) 中给出的误差只是一个粗略的估计。在实际情况中，为了估计达到某一模拟精度所需要的量子电路深度，就需要进一步推导其更严格的误差上界。下面，我们展示一个比较简略的误差上界计算过程，对具体的计算过程不感兴趣的读者可以直接跳至 (11) 阅读这部分的结论。\n",
    "\n",
    "记函数 $f$ 泰勒展开至 $k$ 阶之后的余项为 $\\mathcal{R}_k(f)$, 我们需要用到如下两个结论：\n",
    "\n",
    "$$\n",
    "\\left\\Vert \\mathcal{R}_k \\left( \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right\\Vert\n",
    "\\leq\n",
    "\\mathcal{R}_k \\left( \\exp \\left( \\sum_{k=1}^{L} \\vert \\tau \\vert \\cdot \\Vert h_k \\Vert \\right) \\right),\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\vert \\mathcal{R}_k(\\exp (\\alpha)) \\vert \\leq \\frac{\\vert \\alpha \\vert^{k+1}}{(k+1)!} \\exp ( \\vert \\alpha \\vert ), ~\n",
    "\\forall \\alpha \\in \\mathbb{C}.\n",
    "\\tag{9}\n",
    "$$\n",
    "\n",
    "这两个结论的证明略微有一些繁琐，故在此略去，感兴趣的读者可以参考 [3] 中的 F.1 节。利用 (5) 中的定义，模拟误差可以写为\n",
    "\n",
    "$$\n",
    "\\epsilon\\left(e^{-iH\\tau}, U_{\\rm circuit}\\right) = \\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert.\n",
    "\\tag{10}\n",
    "$$\n",
    "\n",
    "我们已经知道，此时的误差是泰勒展开至一阶后的余项，再利用 (8)，(9) 和三角不等式，可以将 (10) 中的误差上界进行计算：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert\n",
    "=~&\n",
    "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( -i \\sum_{k=1}^{L} h_k \\tau \\right) - \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right \\Vert\n",
    "\\\\\n",
    "\\leq~&\n",
    "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( -i \\sum_{k=1}^{L} h_k \\tau \\right) \\right) \\right \\Vert\n",
    "+\n",
    "\\left \\Vert \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right \\Vert\n",
    "\\\\\n",
    "\\leq~ &\n",
    "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( \\vert \\tau \\vert \\cdot \\left \\Vert \\sum_{k=1}^{L} h_k \\right \\Vert \\right) \\right) \\right \\Vert\n",
    "+ \n",
    "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\sum_{k=1}^{L} \\left( \\vert \\tau \\vert \\cdot \\Vert h_k \\Vert \\right) \\right) \\right \\Vert\n",
    "\\\\\n",
    "\\leq~&\n",
    "2 \\mathcal{R}_1 \\left( \\exp ( \\vert \\tau \\vert L \\Lambda ) \\right )\n",
    "\\\\\n",
    "\\leq~&\n",
    " ( \\vert \\tau \\vert L \\Lambda )^2 \\exp ( \\vert \\tau \\vert L \\Lambda ),\n",
    "\\end{aligned}\n",
    "\\tag{11}\n",
    "$$\n",
    "\n",
    "其中 $\\Lambda = \\max_k \\Vert h_k \\Vert$。随后考虑完整的演化时间 $t = r \\cdot \\tau$，那么模拟长度为 $t$ 的时间演化算符时的误差为：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\left \\Vert \\left ( \\exp\\left(-i\\sum_{k=1}^L h_k \\tau \\right)\\right)^r - \\left (\\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right)^r \\right \\Vert\n",
    "\\leq ~&\n",
    "r \\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert\n",
    "\\\\\n",
    "=~& r (  \\tau L \\Lambda )^2 \\exp ( \\vert \\tau \\vert L \\Lambda )\n",
    "\\\\\n",
    "=~& \\frac{(  t L \\Lambda )^2}{r} \\exp \\left( \\frac{\\vert t \\vert L \\Lambda}{r} \\right).\n",
    "\\end{aligned}\n",
    "\\tag{12}\n",
    "$$\n",
    "\n",
    "这里用到了量子电路中误差线性累积的结论，即 $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$，不熟悉这一结论的读者可以参考 [4] 中的 4.5.3 节。至此，我们就计算出了 product formula 对于一段完整的演化时间 $t$ 的模拟误差上界，即 (7) 式中的二阶项 $O(t^2/r)$。 \n",
    "\n",
    "实际上，我们还可以通过 Suzuki 分解的方法来进一步优化对每一个“时间块”内的时间演化算符 $e^{-iH\\tau}$ 的模拟精度。对于哈密顿量 $H = \\sum_{k=1}^{L} a_k h_k$，其时间演化算符的 Suzuki 分解可以写为\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "S_1(\\tau) &= \\prod_{k=0}^L \\exp ( -i h_k \\tau)，\n",
    "\\\\\n",
    "S_2(\\tau) &= \\prod_{k=0}^L \\exp ( -i h_k \\frac{\\tau}{2})\\prod_{k=L}^0 \\exp ( -i h_k \\frac{\\tau}{2})，\n",
    "\\\\\n",
    "S_{2k}(\\tau) &= [S_{2k - 2}(p_k\\tau)]^2S_{2k - 2}\\left( (1-4p_k)\\tau\\right)[S_{2k - 2}(p_k\\tau)]^2,\n",
    "\\end{aligned}\n",
    "\\tag{13}\n",
    "$$\n",
    "\n",
    "对于 $k > 1, k\\in\\mathbb{Z}$，其中 $p_k = 1/(4 - 4^{1/(2k - 1)})$。先前推导的 product formula 实际上只使用了一阶的 Suzuki 分解 $S_1(\\tau)$ 来对每个“时间块”进行模拟，因此也被称为一阶 Suzuki product formula，或者简称一阶 product formula。Suzuki product formula 在物理中也经常被称为 Trotter-Suzuki 分解方法。对于更高阶的 product formula 而言，使用 (10-12) 中类似的计算，可以证明第 $2k$ 阶 product formula 的整体误差上界为：\n",
    "\n",
    "$$\n",
    "\\epsilon\\left(e^{-iHt}, \\left(S_{2k}(\\tau)\\right)^r\\right)\n",
    "\\leq\n",
    "\\frac{(2L5^{k-1}\\Lambda\\vert t \\vert)^{2k+1}}{3r^{2k}} \\exp \\left( \\frac{2L5^{k-1}\\Lambda\\vert t \\vert}{r} \\right),\n",
    "~ k > 1.\n",
    "\\tag{14}\n",
    "$$\n",
    "\n",
    "在得到了模拟误差上界的基础上，还可以进一步计算达到一定最小精度 $\\epsilon$ 时所需要的电路深度的下界。需要指出的是，(12) 和 (14) 中给出的误差上界的计算是十分宽松的。近年来，许多工作都进一步地给出了更加紧致的上界 [3, 5-6]。此外，也有人提出了不基于 Suzuki 分解的 product formula [7]。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3bf6fb2",
   "metadata": {},
   "source": [
    "![image.png](./figures/trotter_suzuki_circuit.png)\n",
    "<div style=\"text-align:center\">图1：Suzuki product formula 电路的示意图</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0c26758",
   "metadata": {},
   "source": [
    "### 利用 Paddle Quantum 验证基于 Suzuki-product formula 的时间演化电路\n",
    "\n",
    "尽管人们不断的在优化 Suzuki-product formula 的误差上界，但是在实际使用中，真实的误差往往和理论上给出的上界有一定差距，也就是说，我们现在能计算的理论上的 product formula 误差依然只是一个十分宽松的上界 [3]。因此，对于一个真实的物理系统，我们往往需要通过数值实验来计算其真实的误差，从而给出一个经验性的误差上界（empirical bound）。这样的数值试验对于计算在一定精度下模拟某一时间演化过程所需要的电路深度十分重要。\n",
    "\n",
    "在刚刚介绍的 `construct_trotter_circuit` 函数中，用户也可以通过改变传入参数 `tau`、`steps`、`order` 来实现对任意阶 Suzuki product formula 电路的创建。下面，我们就利用该函数的这一功能来验证一下 Suzuki product formula 的误差。\n",
    "\n",
    "继续使用先前的哈密顿量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "86de7992",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H = 1.0 Z0, Z1\n",
      "1.0 X0\n",
      "1.0 X1\n"
     ]
    }
   ],
   "source": [
    "print('H =', h_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e183412",
   "metadata": {},
   "source": [
    "这里我们通过改变 `tau` 和 `steps` 两个参数来将 $t=1$ 的演化过程进行拆分。（提示：$\\tau \\cdot n_{\\rm steps} = t$）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "240ed09f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--*-----------------*----Rx(0.667)----*-----------------*----Rx(0.667)----*-----------------*----Rx(0.667)--\n",
      "  |                 |                 |                 |                 |                 |               \n",
      "--x----Rz(0.667)----x----Rx(0.667)----x----Rz(0.667)----x----Rx(0.667)----x----Rz(0.667)----x----Rx(0.667)--\n",
      "                                                                                                            \n",
      "电路的酉算符和时间演化算符之间的保真度为： 0.984\n"
     ]
    }
   ],
   "source": [
    "# 将长度为 t 的时间演化过程拆为 r 份\n",
    "r = 3\n",
    "t = 1\n",
    "cir = UAnsatz(2)\n",
    "# 搭建时间演化电路，tau 为每个“时间块”的演化时间，故为 t/r，steps 是“时间块”的重复次数 r\n",
    "construct_trotter_circuit(cir, h_2, tau=t/r, steps=r)\n",
    "print(cir)\n",
    "print('电路的酉算符和时间演化算符之间的保真度为： %.3f' \n",
    "      % gate_fidelity(cir.U.numpy(), linalg.expm(-1 * 1j * h_2.construct_h_matrix())))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4708cad9",
   "metadata": {},
   "source": [
    "我们发现，通过将 $t=1$ 的模拟时间拆分为三个“时间块”，模拟误差被成功减小了。\n",
    "\n",
    "如果进一步地增加拆分后的“时间块”的数量的话，误差还可以被进一步减小："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "abc3742d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEJCAYAAACkH0H0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkeklEQVR4nO3deZhcZZn38e9d1XvS6YTsGySQEAhRiSAqDsgrMCwGYdQBouAgjAiOM+IFOODGOCwyig7zjsqIgrhC2EcYEEYY4RUZDEuUVBYSwpJ0ZV+qujvd6e6q+/2jTifd1Us6Sfc5tfw+11VXV51z6py7Kp369XPOU89j7o6IiEhYYlEXICIi5UXBIyIioVLwiIhIqBQ8IiISKgWPiIiEqiLqAorBuHHjfMaMGVGXISJSVF566aUt7j4+f7mCZxBmzJjBiy++GHUZIiJFxcze6mu5TrWJiEioFDwiIhIqBY+IiIRKwSMiIqFS8IiISKjKrlebmY0AfgC0A79z919GXJKISFkpieAxszuBBcAmd5/XbfnpwL8BceDH7n4z8FHgfnd/xMwWAQoeEeHhVxr59hMrSe5oZcroWq4+bQ7nzJ9adjWEUYeVwrQIZnYi0Az8rCt4zCwOvAacCqwDFgMLgbOBx919iZn9yt0/sbf919fX+zHHHNNj2bnnnsvnPvc5du7cyZlnntnrORdddBEXXXQRW7Zs4eMf/3iv9ZdffjnnnXcea9eu5cILL+y1/sorr+Sss85i5cqVfPazn+21/qtf/SqnnHIKS5Ys4Yorrui1/qabbuL444/nD3/4A1/+8pd7rb/11ls5+uij+e1vf8sNN9zQa/0Pf/hD5syZwyOPPMJ3vvOdXut//vOfM336dBYtWsRtt93Wa/3999/PuHHjuOuuu7jrrrt6rX/ssceoq6vjBz/4Affee2+v9b/73e8AuOWWW3j00Ud7rKutreXxxx8H4Prrr+epp57qsX7s2LE88MADAFx77bU8//zzPdZPmzaNX/ziFwBcccUVLFmypMf6ww8/nNtvvx2ASy+9lNdee63H+qOPPppbb70VgAsuuIB169b1WP/+97+fb37zmwB87GMfY+vWrT3Wn3zyyXzta18D4IwzzqC1tbXH+gULFnDVVVcBcNJJJ5FvqH737njqz9z06FIyVaOIt6cZ8/azjNy6IvTfveaxR7D94BN313HVqbO5/Mz3hPq71zz2CLYeejoer9y9XQVZRq/6L0ZuXbF72XD+7jWPPYJts84ga3vaA5bpYOya3+yuIYzfvVknL+TaB1+ltSPTq46LT37nPv3uPfPMMy+5+7H525VEi8fdnzWzGXmLjwNWu/saADO7h1zorAOmAUsY4BqXmV0KXApQXV099EVL2VrdMYYP3Pw0yR2txI/8FKPfeqbHh1sYHn6lkW/9zzoy1Q0AZKob2Hro6UO2f3cn645bHLcYmOHEwGJsa82Q3NHK5lZnx6RjSB18Ah6r3F3Hd5/bzK6q14iloW3k1NxzzQADM154q4nVLZtYnoqzc8wsPFgOue0eX7aFESPbeGVHNU3j54HFdq8D4+cvrCNeUcEf0/WkJh+HG6SmvLdH6AB0EmPrzL+kvW4CGIBRUVHJ9Y8uI+vOC7umsXXGyQQrcTOaa2q45oE/k3Vnsc9m82Gjd6/HjJa6Oj7/q5dxh1cq38mO2YcE7w20NczAredHsscr2XLoGTRPeCcAv6+s5/zbn8cdlo8/mV2j2nfvG+BJr2f595/DgVXTz6JzcqbH+oda63nu1mdxhzcPX0i2q+ER1PDzbSNpvu9PZLI9GyQer2T7wScCO/bn16GXkmjxAATB82i3Fs/HgdPd/W+DxxcC7wX+Efge0Ab8fjDXeI499ljXyAUyFB5+pbHXX5O1lXG++dF37D6Vkc06HdksHRmnvTNLRyZLe2eW9sye+7mfnlvW9Xj3Oqe9M5P72WP7ru2ch19p7FFDl+qKGEdPH00m63Rmc+HRmfHgcZasQ2c2SybTbX3WyWScTNf94FYqquIxLMi2mBkxM4zgcSx3P2YWbGPEDIzgZ7A8lvez6zndn7d8fbrfGo45ZMzuYxq5HfQ8bm55kC89jmPd6u1xn97PNTMe+VOyzxoMeOPmD+/Te2dmpdvi2Rfu3gJ8Ouo6JHxDdd46m3V2dmRo2dVJy65OdrZn9vxs7wyWZ9jZ3klLe4adu3I/W3Z18vSKTezqzPbYX2tHhi/eu4QvP/QqHZlccAy1yrhRGY9RGY9RVRHrM3SA3bVVV8aoNaMiZsRjseBn7lYRM2Ix272s5+NYr/XxAba96r4/9VvzTy8+jli3D/uYQTxmuz/c4zHbsy4GcetjXSx4nLfOzIJtcvv/0C2/I5lq61XD1NG1PHfNh4bmH2EvPnDz0zTuaO21fOroWh64/PhQagB4+a3tfdYxZXTtkB2jlIOnEZje7fG0YJmUoYdfaeSaB/9MW0fug7VxRytX3/8nnl+zhTkTR7GzvZPmrrDIC43mIFS61vX3od2XyrgxorqCEVUV1FXFe4VOF3f4xHEHU1WxJxyq4jEq40ZVRTz42bUsRmVwv6rCdm9fGe9aFuu2zKiMxYjFrMfxBvqQW/TZ9+/DO3tg/vW/X+u3jg8e3mtsyWHzpdOP6LMlevVpc0Kr4erT5kReQ1h1lHLwLAZmm9lMcoFzPrDXjgRSnNo6MmxItZFMtbIh1cb6VBvrU62s35G7v2JDmvyzPx0ZZ9HiPRdnq+Ix6qrju0NiRHUFI6rjjBlRx4iqOHXVFYzoWl5V0ce2FT22q6uqoKqi52XEgT7wv7pg7rC8N30ppw+5wehq+UbZo6wQagirjpK4xmNmdwMnAeOAjcB17n6HmZ0J3EquO/Wd7n7j/uxf13gO3IGc5uovVDak2kjuaGNDuo1tLe29njemrpJJDbVMbqjh6RWb+ty3AUu+/pfUVsV7hcRwGMw1nrCUS9ddiU5/13hKIniGm4LnwAz0YXv6vEk9WyiDDJXRdZVMDkJlzy14PLqWSaNqqK2K795+oJZGWOfwu+iDVsqFgucAKHgOTH8f+jGj1+kvyIXKpFE1TBldy6SGGqY01DCpoTb4mQuY7qEyGIXU0hApF+rVJqFr78zy1PKNfYYO5ELnylMPZ/LoPS2XSQ011FUN/a9loZw/FxEFjwyD1ZuauffFtTzw0jq2trT327KZOrqWvz95dmh1nTN/qoJGpAAoeGRItLZn+K9X17No8dssfnM7FTHjlCMnct5x09ne3M5XHl4aec8lESkMCh45IEsbU9yz+G3+85UkTbs6mTluBNeccQQfffdUJtTX7N4uFjOd5hIRQMEj+yHV2sGvlzRyz+K1JJJpqitifPgdkznvPdM5buZBmFmv5+g0l4h0UfDIoLg7f3xjG4sWr+W/Xl3Prs4scyeP4vqzj+IjR0+lobZy7zsREUHBI3uxuWkXD768jkWL17JmSwv11RV8/JhpLDzuYOZNbYi6PBEpQgoe6SWTdZ5dtZlFf1zLb5dvpDPrvGfGGD73f2Zx5jsmDUt3ZxEpH/oEkd3Wbd/JvS+u474X17I+1cbYEVVc/BczOffY6cyaMDLq8kSkRCh4ylx7Z5b/XraRexa/ze9XbwHghNnj+dqCuZxy5MRQxi8TkfKi4ClTqzc1sWjxWh54uZFtLe1MaajhHz40m78+dhrTxtRFXZ6IlDAFT4nrPiDlpIYaTpoznlUbm3nxrdyXPE+dO5Hz3jOdE2aPJx7r3Q1aRGSoKXhKWP7AmOtTbdz9x7WMr6/iy2cewV/Nn8b4+uqIqxSRcqPgKWHffmJln7NlVsVjXHriYRFUJCICunJcwpL9jAqd3NF7bnkRkbAoeErYlNG1+7RcRCQMCp4SdvVpc3p1GNCo0CISNQVPCTtn/lSmj6mlKh7DyM1/oxk3RSRq6lxQwjJZZ1PTLhYeN51vnD0v6nJERAC1eEram1tb2Nme4agpGsxTRAqHgqeEJZJpAI6aOiriSkRE9lDwlLBEMkVl3Jg9oT7qUkREdlPwlLBlyTSHT6zXQJ8iUlD0iVSi3J1EMs1RU3SaTUQKi4KnRG1It7GtpV0dC0Sk4Ch4StTSxqBjgVo8IlJgFDwlKpFMYQZHTlbwiEhhUfCUqEQyzcyxIxhRre8Ii0hhUfCUqGXJNHN1mk1ECpCCpwRtb2mncUerOhaISEFS8JSgZevVsUBECpeCpwQlkilAwSMihUnBU4ISyTSTRtUwdmR11KWIiPSi4ClBiWSaeRoYVEQKlIKnxLS2Z1izuZm56lggIgVKwVNilm9Ik3Vd3xGRwqXgKTG75+BR8IhIgVLwlJhlyRQNtZVMHV0bdSkiIn1S8JSYrqkQzCzqUkRE+qTgKSEdmSwr1jfpNJuIFDQFTwlZvamZ9kxWQ+WISEFT8JQQdSwQkWKg4CkhiWSKmsoYh44fGXUpIiL9KtvgMbNDzewOM7s/6lqGSiKZ5ohJo4jH1LFARApXaMFjZl8ws6VmljCzKw5gP3ea2SYzW9rHutPNbKWZrTazawbaj7uvcfdL9reOQpPNOsuDHm0iIoUslOAxs3nAZ4DjgHcBC8xsVt42E8ysPm9Zj20CdwGn93GMOPB94AxgLrDQzOaa2TvM7NG824QheWEFZO32nTTt6lTHAhEpeGG1eI4EXnD3ne7eCTwDfDRvmw8CD5tZNYCZfQb49/wdufuzwLY+jnEcsDpoybQD9wBnu/ur7r4g77ZpMEWb2VlmdnsqlRr0C41KV8cCDQ4qIoUurOBZCpxgZmPNrA44E5jefQN3vw94AlhkZp8ELgb+eh+OMRVY2+3xumBZn4Ja/gOYb2bX9rWNuz/i7pc2NBR+KyKRTBGPGYdPrN/7xiIiEaoI4yDuvtzM/gV4EmgBlgCZPrb7lpndA9wGHObuzcNY01bgsuHaf9gSyTSzJ4ykpjIedSkiIgMKrXOBu9/h7se4+4nAduC1/G3M7ARgHvAQcN0+HqKRnq2oacGyspBIppmrjgUiUgTC7NU2Ifh5MLnrO7/KWz8fuB04G/g0MNbMbtiHQywGZpvZTDOrAs4Hfj0UtRe6TU1tbG7apY4FIlIUwvwezwNmtgx4BPg7d9+Rt74OONfdX3f3LPAp4K38nZjZ3cDzwBwzW2dmlwAEnRY+T+460XLgXndPDNurKSAasUBEikko13gA3P2Evax/Lu9xB/CjPrZbOMA+HgMe298ai9WyIHh0qk1EikHZjlxQSpY2pjj4oDpG1VRGXYqIyF4peEpAQiMWiEgRUfAUuXRbB29v26ngEZGioeApcst2dyxQjzYRKQ4KniKnHm0iUmwUPEUukUwxbmQ1E0bVRF2KiMigKHiK3LJkWgODikhRUfAUsbaODKs2Nes0m4gUFQVPEXttYxOZrKtjgYgUFQVPEVPHAhEpRgqeIpZIpqivrmD6mLqoSxERGTQFTxFLJNMcOWUUsZhFXYqIyKApeIpUJuusWN+k02wiUnQUPEXqjS3NtHZk1LFARIqOgqdILW1UxwIRKU4KniKVSKaoqogxa8LIqEsREdknCp4ilUimmTOxnsq4/glFpLjoU6sIubvm4BGRoqXgKUKNO1pJtXYoeESkKCl4itDuEQumqkebiBQfBU8RSiTTxAyOnKQWj4gUn0EFj5k9ZGbnmFnlcBcke7csmeLQ8SOprYpHXYqIyD4bbIvn/wFfBzaY2W1mdvww1iR7oY4FIlLMBhU87v5dd383cCKwA7jbzFaZ2dfN7LDhLFB62tbSzvpUm4JHRIrWPl3jcfeEu18LXADsBK4DXjaz35rZu4ajQOkpkUwBaKgcESlagw4eM5tjZteb2evA7cAiYAYwEXgMeHg4CpSeNAePiBS7isFsZGYvkguZRcAn3P2FvE2+a2Z/P8S1SR8SyTRTR9cyuq4q6lJERPbLoIIHuBn4tbu397eBu88cmpJkIIlkirlq7YhIERvsqbav9BU6QUtIQtKyq5M3trToNJuIFLXBBk+vnmtmZsChQ1uODGT5+jTu6lggIsVtwFNtZvaz4G51t/tdZgCJ4ShK+qaOBSJSCvZ2jef1fu478Bxw35BXJP1KJFOMqatkckNN1KWIiOy3AYPH3b8BYGb/6+5PhFOS9CeRTDNvagO5s5wiIsWp3+AxsxPd/dngYYeZfaiv7dz96WGpTHpo78zy2sYmLv4LdR4UkeI2UIvnB8C84P4d/WzjqINBKFZtaqIj4+pYICJFr9/gcfd53e7rz+yIqWOBiJQKzcdTJJYl09RVxZk5dkTUpYiIHJCBrvGsJXcqbUDufvCQViR9SiRTHDl5FLGYOhaISHEb6BrPBaFVIQPKZp1lyTQfO2Za1KWIiBywga7xPBNmIdK/t7btpKU9o+s7IlISBjv1dbWZ3Whma8wsFSz7SzP7/PCWJ6A5eESktAy2c8G/kuta/Un2XPdJAJcPR1HSUyKZpiJmzJ44MupSREQO2GCnRfgrYJa7t5hZFsDdG81s6vCVJl2WNqaYPbGe6op41KWIiBywwbZ42skLKTMbD2wd8oqkB/dcxwJd3xGRUjHY4LkP+KmZzQQws8nA94B7hqswydmY3sXWlnYFj4iUjMEGz5eBN4BXgdHAKiAJfGN4ypIuXR0L5k1VxwIRKQ2DusYTzD76ReCLwSm2Le6+1y+XyoFLJNOYwZGT1eIRkdIw0MgFAw3+Wd81NL+7rxnqomSPRDLFjLEjGFk92H4gIiKFbaBPs9Xkuk4be7pQd43X0r21o65WwyiRTPOu6aOjLkNEZMj0e43H3WPuHnf3GPC35DoSzAFqgCOAXwGXhFJlmUrt7GDd9lZ1LBCRkjLY8zfXA7PdvTV4vMrMPgu8Btw1HIUJJNZrxAIRKT2D7dUWA2bkLTuEIj7NZmaHmtkdZnZ/1LX0Z5nm4BGRErQvQ+Y8bWY3mdnlZnYT8FSwfFDM7ItmljCzpWZ2t5nV7E/BZnanmW0ys6V9rDvdzFaa2Wozu2ag/bj7Gncv6FOFiWSaiaOqGTeyOupSRESGzKCCx92/DXwamAh8BJgEXOzu3xrM84Ohdf4BODaY2TQOnJ+3zQQzq89bNquP3d0FnN7HMeLA94EzgLnAQjOba2bvMLNH824TBlN31BLJlE6ziUjJGXQfXXf/DfCbAzxWrZl1AHXkvoDa3QeBy8zsTHffZWafAT5KLki61/Gsmc3oY//HAau7uneb2T3A2e7+TWDB/hRsZmcBZ82a1Vf+Da+2jgyvb27htKMmhX5sEZHhNND3eL7i7jcG9/+5v+3c/et7O0gwoOgtwNtAK/Ckuz+Zt819wZA8i8zsPuBi4NTBvQwApgJruz1eB7y3v43NbCxwIzDfzK4NAiq/7keAR4499tjP7EMdQ2LFhiYyWdf1HREpOQO1eL5B7oMZ4DByA4XuFzMbA5wNzAR2APeZ2QXu/ovu27n7t4KWym3AYe7evL/H3Bt33wpcNlz7P1BLG9WjTURK00DBs7Pb/bPc/UD+9D4FeMPdNwOY2YPA8UCP4DGzE8jN+/MQcB2wLxPNNQLTuz2eFiwrSolkmlE1FUwbUxt1KSIiQ2qg4HndzL5DbsK3CjP7NHtGLtjN3e8cxHHeBt5nZnXkTrWdDLzYfQMzmw/cTu56zBvAL83sBnf/6qBeCSwGZgen6xrJdV74xCCfW3CWBR0LuoYmEhEpFQMFz3nAl4CFQBXwqT62cWCvwePuLwTfl3kZ6AReIRcy3dUB57r76wBm9ingovx9mdndwEnAODNbB1zn7ne4e2cwFfcT5HrN3enuib3VVog6M1lWbGjiwvcdEnUpIiJDrt/gcffXyA2Vg5k95e4nH8iB3P06cqfP+lv/XN7jDuBHfWy3cIB9PAY8dgBlFoTXN7ewqzPLUVPVsUBESs9gv8dzQKEj+6ZrDh51LBCRUjTYkQskRIlkmuqKGIeOGxF1KSIiQ07BU4ASyRRHTB5FRVz/PCJSevTJVmDcnWXJtL44KiIlS8FTYNZtbyXd1qngEZGSpeApMOpYICKlTsFTYBLJNPGYccSk+r1vLCJShBQ8BSaRTHPY+BHUVBbtHHsiIgNS8BSYpY2ag0dESpuCp4BsbtrFpqZd6lggIiVNwVNA1LFARMqBgqeAJJJpAOaqxSMiJUzBU0CWJdNMP6iWhtrKqEsRERk2Cp4CkkimOGqyTrOJSGlT8BSIprYO3ty6Ux0LRKTkKXgKxPL1TQCag0dESp6Cp0CoR5uIlAsFT4FIJNOMG1nFhPrqqEsRERlWCp4CkUimmTulATOLuhQRkWGl4CkAuzozrNrYpI4FIlIWFDwFYNXGZjqzruARkbKg4CkA6lggIuVEwVMAljamGVldwSEH1UVdiojIsFPwFIBEMsXcyaOIxdSxQERKn4InYpmss3x9kwYGFZGyoeCJ2BtbWmjtyKhjgYiUDQVPxNSxQETKjYInYsuSaariMWZPHBl1KSIioVDwRCyRTHP4pJFUxvVPISLlQZ92EXJ3zcEjImVHwROh9ak2tu/s0FQIIlJWFDwRSiTTAOrRJiJlRcEToUQyhRkcMUnBIyLlQ8EToUQyzcxxIxhRXRF1KSIioVHwRGhZMq3v74hI2VHwRGR7SzuNO1p1fUdEyo6CJyJdHQvmqcUjImVGwRORPUPlqMUjIuVFwRORRDLNlIYaxoyoiroUEZFQKXgikkimmKvTbCJShhQ8EdjZ3smaLS06zSYiZUnBE4Hl65tw1/UdESlPCp4ILOvqWDBVp9pEpPwoeCKQSKYZXVfJlIaaqEsREQmdgicCiWSao6aMwsyiLkVEJHQKnpB1ZLKs3NCkoXJEpGwpeEK2elMz7ZmsOhaISNlS8IRMc/CISLlT8IQskUxRWxln5riRUZciIhIJBU/IEsk0R06uJx5TxwIRKU8KnhBls645eESk7Cl4QvT2tp007+rU9R0RKWsKnhDt6VigFo+IlC8FT4gSyRQVMePwSepYICLlS8ETokQyzawJI6muiEddiohIZBQ8IUqoY4GIiIInLJvSbWxp3qWOBSJS9hQ8IdGIBSIiOWUbPGZ2qJndYWb3h3G8RDAHz1wFj4iUuVCCx8zmmNmSbre0mV2xn/u608w2mdnSPtadbmYrzWy1mV0z0H7cfY27X7I/NeyPRDLNIWPrqK+pDOuQIiIFqSKMg7j7SuBoADOLA43AQ923MbMJQKu7N3VbNsvdV+ft7i7ge8DP8p4fB74PnAqsAxab2a+BOPDNvH1c7O6bDuxV7ZtEMs28qWrtiIhEcartZOB1d38rb/kHgYfNrBrAzD4D/Hv+k939WWBbH/s9DlgdtGTagXuAs939VXdfkHcLNXTSbR28vW2nerSJiBBN8JwP3J2/0N3vA54AFpnZJ4GLgb/eh/1OBdZ2e7wuWNYnMxtrZv8BzDeza/vZ5iwzuz2VSu1DGb0tU8cCEZHdQg0eM6sCPgLc19d6d/8W0AbcBnzE3ZuHqxZ33+rul7n7Ye6efyqua5tH3P3ShoYDa6loqBwRkT3CbvGcAbzs7hv7WmlmJwDzyF3/uW4f990ITO/2eFqwLHKJxhQT6qsZX18ddSkiIpELO3gW0sdpNgAzmw/cDpwNfBoYa2Y37MO+FwOzzWxm0LI6H/j1AdY7JHIjFug0m4gIhBg8ZjaCXI+zB/vZpA44191fd/cs8CkgvwMCZnY38Dwwx8zWmdklAO7eCXye3HWi5cC97p4Y+leyb9o6Mqze3KzTbCIigVC6UwO4ewswdoD1z+U97gB+1Md2CwfYx2PAYwdQ5pBbuaGJTNbV4hERCZTtyAVhUccCEZGeFDzDLJFMUV9TwfSDaqMuRUSkICh4hlkimWbu5FGYWdSliIgUBAXPMMpknRUbNAePiEh3Cp5htGZzM20dWXUsEBHpRsEzTB5+pZFzf/g8AP/ymxU8/EpBfJdVRCRyoXWnLicPv9LItQ++SmtHBoBNTbu49sFXAThnfr/Dx4mIlAW1eIbBt59YuTt0urR2ZPj2EysjqkhEpHAoeIZBckfrPi0XESknCp5hMGV039/Z6W+5iEg5UfAMg6tPm0NtZbzHstrKOFefNieiikRECoc6FwyDrg4E335iJckdrUwZXcvVp81RxwIRERQ8w+ac+VMVNCIifdCpNhERCZWCR0REQqXgERGRUCl4REQkVAoeEREJlbl71DUUPDPbDLwVdR0HaBywJeoiCoTei570fvSk92OPA30vDnH38fkLFTxlwsxedPdjo66jEOi96EnvR096P/YYrvdCp9pERCRUCh4REQmVgqd83B51AQVE70VPej960vuxx7C8F7rGIyIioVKLR0REQqXgERGRUCl4SpiZTTez/zGzZWaWMLMvRF1TITCzuJm9YmaPRl1L1MxstJndb2YrzGy5mb0/6pqiYmZfDP6fLDWzu82sJuqawmRmd5rZJjNb2m3ZQWb232a2Kvg5ZiiOpeApbZ3Ale4+F3gf8HdmNjfimgrBF4DlURdRIP4N+I27HwG8izJ9X8xsKvAPwLHuPg+IA+dHW1Xo7gJOz1t2DfCUu88GngoeHzAFTwlz9/Xu/nJwv4nch0pZTxJkZtOADwM/jrqWqJlZA3AicAeAu7e7+45Ii4pWBVBrZhVAHZCMuJ5QufuzwLa8xWcDPw3u/xQ4ZyiOpeApE2Y2A5gPvBBxKVG7FfgSkI24jkIwE9gM/CQ49fhjMxsRdVFRcPdG4BbgbWA9kHL3J6OtqiBMdPf1wf0NwMSh2KmCpwyY2UjgAeAKd09HXU9UzGwBsMndX4q6lgJRAbwbuM3d5wMtDNGplGITXLs4m1wYTwFGmNkF0VZVWDz33Zsh+f6NgqfEmVkludD5pbs/GHU9EfsA8BEzexO4B/iQmf0i2pIitQ5Y5+5dreD7yQVROToFeMPdN7t7B/AgcHzENRWCjWY2GSD4uWkodqrgKWFmZuTO3y939+9GXU/U3P1ad5/m7jPIXTh+2t3L9q9ad98ArDWzOcGik4FlEZYUpbeB95lZXfD/5mTKtKNFnl8DfxPc/xvgP4dipwqe0vYB4EJyf9kvCW5nRl2UFJS/B35pZn8GjgZuiracaAStvvuBl4FXyX02ltXQOWZ2N/A8MMfM1pnZJcDNwKlmtopcq/DmITmWhswREZEwqcUjIiKhUvCIiEioFDwiIhIqBY+IiIRKwSMiIqFS8IiExMzeNLNTIjr2RDN71syazOw7UdQg0qUi6gJEJBSXAluAUb4P36Ews5OAX7j7tGGqS8qQWjwiRSYYPXlfHQIs25fQERkuCh4pa8Hpr6vM7M9mljKzRV0TgJnZRWb2+7zt3cxmBffvMrMfmNnjZtZsZs+Z2SQzu9XMtgeTq83PO+R7gon5tpvZT7pPNmZmC4LRJXaY2R/M7J15df5jMMJAS1/hY2bHm9ni4HUsNrPju+okN9zJl4I6e53uM7Mzg7qazKwxeE9GAI8DU4LnNZvZFDOLmdk1Zva6mW01s3vN7KBgPzOC9+hSM0ua2Xozu6rbcY4zsxfNLG1mG82s7IdyKkvurptuZXsD3gT+SG5E4oPIjc91WbDuIuD3eds7MCu4fxe501fHADXA08AbwKfITSR2A/A/ecdaCkwPjvUccEOwbj65ARjfGzz3b4Ltq7s9d0nw3No+XsdBwHZyQyRVAAuDx2O71XrDAO/DeuCE4P4Y4N3B/ZPIDSTafdsvAP8LTAOqgR8CdwfrZgTv0d3ACOAd5KZeOCVY/zxwYXB/JPC+qH8HdAv/phaPCPxfd0+6+zbgEXJjlg3WQ+7+kru3AQ8Bbe7+M3fPAIvIBUp333P3tcGxbiQXEJC7BvNDd3/B3TPu/lNgF7mZY7vXudbdW/uo48PAKnf/ubt3uvvdwArgrEG+jg5grpmNcvftHkwg2I/LgK+4+zp33wX8E/DxvFbYN9y9xd1fBX7S7XV2ALPMbJy7N7v7/w6yPikhCh6R3ARXXXaS+0t8sDZ2u9/ax+P8fa3tdv8tci0tyF2DuTI4zbbDzHaQa91M6ee5+aYE++vuLQY/4+zHgDOBt8zsGTN7/wDbHgI81K3O5UCGnpOE9fc6LwEOB1YEpwMXDLI+KSEKHpH+tZCbAhkAM5s0BPuc3u3+weyZXnktcKO7j+52qwtaLl0G6hiQJBcI3R0MNA6mKHdf7O5nAxOAh4F7BzjmWuCMvFprPDeLZ5c+X6e7r3L3hcFx/gW4v1xnPS1nCh6R/v0JOMrMjg46AfzTEOzz78xsWnAx/ivkTscB/Ai4zMzeazkjzOzDZlY/yP0+BhxuZp8wswozOw+YCzy6tyeaWZWZfdLMGjw3CVqaPVODbwTGmllDt6f8B3CjmR0SPH+8mZ2dt9uvBXPbHAV8uut1mtkFZjbe3bPAjmBbTUNeZhQ8Iv1w99eAfwZ+C6wCfj/wMwblV8CTwBrgdXIdEHD3F4HPAN8j1ylgNbnODYOtdSuwALgS2Ap8CVjg7lsGuYsLgTfNLE3uGs4ng/2uINdRYE1wam0K8G/kJgh70syayHU0eG/e/p4JXsNTwC3u/mSw/HQgYWbNwX7O7+ealZQwzccjIkPGzGaQ69lX6e6dEZcjBUotHhERCZWCR0REQqVTbSIiEiq1eEREJFQKHhERCZWCR0REQqXgERGRUCl4REQkVP8f1SBySEXruNgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 导入一些作图和计算误差需要的额外包\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def get_fid(n_steps):\n",
    "    t = 1\n",
    "    cir = UAnsatz(2)\n",
    "    construct_trotter_circuit(cir, h_2, tau=t/n_steps, steps=n_steps)\n",
    "    return gate_fidelity(cir.U.numpy(), linalg.expm(-1 * 1j * h_2.construct_h_matrix()))\n",
    "plt.axhline(1, ls='--', color='black')\n",
    "plt.semilogy(np.arange(1, 11), [get_fid(r) for r in np.arange(1, 11)], 'o-')\n",
    "plt.xlabel('number of steps', fontsize=12)\n",
    "plt.ylabel('fidelity', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8381126",
   "metadata": {},
   "source": [
    "不仅如此，我们还可以通过改变 product formula 的阶数来减小模拟误差。目前，`construct_trotter_circuit()` 函数支持通过参数 `order` 来实现任意阶数的 Suzuki product formula。下面，就让我们来分别计算一下一阶和二阶时间演化电路的误差随着 $\\tau$ 大小的变化，并和上文中计算的理论误差上界进行对比："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d3256b18",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAENCAYAAADOhVhvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyyElEQVR4nO3dd3zUVb7/8ddJSEKvoSUBCWJClQQDqKiADSw0O3ZxLevqXXdXr7rXe3WLV/fq+nPdZVWsawMVMYCromIBGzWhCUiHTCBAIKGlz/n98Z1oiJTJZCbfmcz7+XjkAfnOzDcfwpAP55zP+RxjrUVERCQQMW4HICIikUtJREREAqYkIiIiAVMSERGRgCmJiIhIwJREREQkYE3cDqChJSYm2h49ergdhohIxFiyZMlua23HIz0WNUnEGDMGGNOrVy8WL17sdjgiIhHDGLPlaI9FzXSWtXa2tfbWNm3auB2KiEijETVJREREgi9qkogxZowxZkpxcbHboYiINBpRsyZirZ0NzM7Kyrql9mMVFRXk5eVRWlrqQmTho2nTpqSkpBAXF+d2KCISIaImidRcWK8tLy+PVq1a0aNHD4wxDR9cGLDWUlhYSF5eHqmpqW6HIyIRImqms461sF5aWkqHDh2iNoEAGGPo0KFD1I/GRKRuoiaJHE80J5Bq+h6INFIle+H7WSG5ddRMZ4mIRKWKUph6NXiWQPJSaJMS1NtHTRI51ppIXWXneHh8zlryi0pIatuMe0elMz4zuf5BHkNVVRWxsbFH/fxIrLVYa4mJ0YBTJCp5q2DGL2DrN3Dpi0FPIBBF01nB2myYnePhgRkr8BSVYAFPUQkPzFhBdo6nXvd9/fXXGTJkCBkZGdx2221UVVXRsmVLfve73zFw4EC+/fbbn33+5JNP0r9/f/r3789TTz0FwObNm0lPT+f666+nf//+bNu2rV5xiUiEshY+vA9Wz4ZRj8KAy0LyZaJmJOKvP8xexff5+476eM7WIsqrvIddK6mo4j+nL2fqwq1HfE3fpNY8NKbfUe+5evVq3nrrLb7++mvi4uK44447eOONNzh48CBDhw7lr3/9K8Bhny9ZsoSXX36ZBQsWYK1l6NChDB8+nHbt2rFu3Tr+9a9/ceqppwbwHRCRRmH+X2HR83D6XXDaHSH7MlGTRII1nVU7gRzvuj/mzp3LkiVLGDx4MAAlJSV06tSJ2NhYLr300h+fV/Pzr776igkTJtCiRQsALrnkEubPn8/YsWM54YQTlEBEolnOG/DZn2DAFXDuH0P6paImiRxrs2FNxxoxAAx77DM8RSU/u57cthlv3XZaoLFxww038Oijjx52/Yknnjhs3aNp06bHXQcBfkwsIhKF1n0Cs+6CniNg3GQI8Zpo1KyJBMu9o9JpFnf4D/JmcbHcOyo94Huec845TJ8+nZ07dwKwZ88etmw5atNMAM4880yys7M5dOgQBw8e5L333uPMM88MOAYRaQQ8S+Dt66FzX7jiNWgSH/IvGTUjkWCprsIKZnVW3759+fOf/8z555+P1+slLi6OyZMnH/M1gwYN4sYbb2TIkCEA/OIXvyAzM5PNmzcHHIeIRLDCDfDGFdCiI1zzLjRt3SBf1lhrG+QLhYusrCxb+zyR1atX06dPH5ciCi/6XohEoAM74YVzofwATPoYEuu/laEmY8wSa23WkR6LmuksdfEVkUapbD+8cRkc3AVXvx30BHI8UZNEdCiViDQ6leXOGsiOlXD5K5ByxMFCSGlNREQkElnrVGFt+AzG/gPSRrkSRtSMREREGpVPH4bl02DkgzDoOtfCUBIREYk03z0LXz8FWZPgrHtcDUVJREQkkqycAR/dD70vhgufAJePcIiaJBLu1VmbN2+mf//+Dfb1Hn74YZ544okG+3oiEgSb5sN7t0G3oXDpCxBz/A4WoRY1SUTVWSIS0QpWwbRroF0qTJwKcc3cjgiIoiQSCSorK7nmmmvo06cPl112GYcOHWLu3LlkZmYyYMAAJk2aRFlZGQA9evRg9+7dACxevJgRI0YAzghj0qRJjBgxgp49e/L000//eP9HHnmEtLQ0zjjjDNauXdvgfz4RCVDRNnj9UohvAde+C83bux3Rj1TiW9uH98OOFcG9Z5cBcMFjx33a2rVrefHFFxk2bBiTJk3iySef5LnnnmPu3LmkpaVx/fXX88wzz3D33Xcf8z5r1qzh888/Z//+/aSnp/PLX/6S5cuXM23aNHJzc6msrGTQoEGccsopQfoDikjIHNrjJJDyQzDpQ2jbze2IDqORSBjp1q0bw4YNA+Daa69l7ty5pKamkpaWBsANN9zAvHnzjnufiy66iISEBBITE+nUqRMFBQXMnz+fCRMm0Lx5c1q3bs3YsWND+mcRkSCoKIGpV8HeTXDVG9D52F3G3aCRSG1+jBhCxdSqsmjbti2FhYVHfG6TJk3wep0zTEpLSw97LCEh4cffx8bGUllZGeRIRSTkvFXw7i9g20K4/GVIDc8u3RqJhJGtW7fy7bffAvDmm2+SlZXF5s2bWb9+PQCvvfYaw4cPB5w1kSVLlgDw7rvvHvfeZ511FtnZ2ZSUlLB//35mz54doj+FiNSbtfDBPbDmfbjgL9BvgtsRHZWSSBhJT09n8uTJ9OnTh7179/Kb3/yGl19+mcsvv5wBAwYQExPD7bffDsBDDz3Er3/9a7Kysvw6qGrQoEFceeWVDBw4kAsuuODHUxRFJAzNewIWvwTD7oaht7kdzTFFTSv4Gsfj3rJu3brDHlP785/oeyHisqWvOj2xTr4KJjzr+mZCUCt4QPtERCQC/DAHZt8NJ54N4/4RFgnkeKImiYiIhLW8xfD2Dc6WgCtehdg4tyPyi5KIiIjbdq+HNy6HVp3hmncgoZXbEflNScQnWtaGjkXfAxEX7C+A1yeAiYFrZ0DLTm5HVCdKIkDTpk0pLCyM6h+i1loKCwtp2rSp26GIRI/SffDGpXBwN1zzNnQ40e2I6kybDYGUlBTy8vLYtWuX26G4qmnTpqSkpLgdhkh0qCyHt66Fgu/h6rcgOTLbECmJAHFxcaSmprodhohEC68XZt4Bm76E8c/ASee5HVHANJ0lItLQPv0fWPEOnPM/kHG129HUi5KIiEhD+nYyfPN3GHwLnPFbt6Opt4hOIsaYnsaYF40x092ORUTkuFZMhzm/hz5jnZ5YEbCZ8HjCLokYY14yxuw0xqysdX20MWatMWa9MeZ+AGvtRmvtze5EKiJSBxu/hPduh+6nwyXPh8XRtsEQdkkEeAUYXfOCMSYWmAxcAPQFJhpj+jZ8aCIiAdixwjnatkMvmPgmxDWeUvqwSyLW2nnAnlqXhwDrfSOPcmAaMK7BgxMRqau9W+D1y6Bpa+do22bt3I4oqMIuiRxFMrCtxud5QLIxpoMx5lkg0xjzwNFebIy51Riz2BizONr3gohIA6o+2rayxEkgbZLdjijoInqfiLW2ELjdj+dNAaYAZGVlRe+2dBFpOOWH4M0roGgrXJ8NnRrnEQuRMhLxADVPp0/xXfObMWaMMWZKcXFxUAMTEfmZqkqYPsnpzHvpC3DC6W5HFDKRkkQWAScZY1KNMfHAVcCsutxA54mISIOwFv79W/jhQ7jwceg71u2IQirskogxZirwLZBujMkzxtxsra0E7gTmAKuBt621q+p4X41ERCT0vvwLLP0XnPk7GHKL29GEXNQcj1stKyvLLl682O0wRKQxWvwyvH83ZFwD4yY3is2EoONxRURCb80HzjRWr/NgzN8aTQI5nqhJIprOEpGQWTnDWUjvmgGXvxIxR9sGQ9QkES2si0jQlRbDjFth+k3QuS9c/TYktHQ7qgYV0ftERERcs/krpxfWvnwYfj+cdU9UjUCqRU0SMcaMAcb06tXL7VBEJJJVlsHnj8DXT0P7VLj5Y0g54ppzVNB0loiIv3auhufPga//BqfcALfNj+oEAlE0EhERCZjXCwuehU8fhoRWMHEapF/gdlRh4bhJxBhzfR3vmWutXR5gPCIi4aXYA9m/dM5DTxsNY/8OLTu5HVXY8GckMhLwZ0ei8T2vCAi7JKI1ERGps5Xvwvu/gaoKuPgpOOXGqNn/4S9/kshmnOTgz3euOomEHWvtbGB2VlZW4+9DICL1U1IEH9wLK96G5Cy4ZAp0ONHtqMKSv2si/qZepWgRiWyb5julu/u3w4gH4Mx7IFbLx0dz3O+MtfYPDRGIiIirKsvgsz/DN39X6W4d+J1ejTGdgFHAQKAtzrTVMuATa+2OUAQXTFoTEZGjKvje2XlesAJOuQlGPQLxLdyOKiIcd5+IMaaPMWY6Tgv264A4YIfv1+uAVcaY6caYviGNtJ60T0REfsbrhW8nw5QRcGCHU7o75iklkDrwZyTyCvA4cI21tqz2g8aYBGAs8CJwWlCjExEJlcNKdy/wle52dDuqiOPPmsjQ6t8bY2Kstd5aj5cB7/g+RETC34rpTtv2qkqnbfugG1S6G6C6lhwUAErVIhKZVLobdH4lEWPMQOB7oOlRHt9qre0ezMBERIJKpbsh4e938H2gExDjOwM9F6cyKxdncT7sV6tVnSUSpSrL4LM/wTf/gPY9VbobZH518bXWdgOSgQpgPtAT+AOwAdgGvB6qAINF1VkiUahgFTx/trP345Qb4XZ13Q02v8dy1trdxpgB1toN1deMMQZoZq09FJLoREQC4fXCd/+EuX+Apm1g4luQPtrtqBolf9dEnsOZusoxxmyvThrWWgsogYhI+CjO85XuzoP0C2HM0yrdDSF/RyKlwOXAw0B7Y8x6fEnF92uutXZnCOITEfHfYaW7T8Og61W6G2J+JRFr7a+rf2+M6Yyze30jkAX8EugOxIYiQBGR4yopgg/ugRXvQMpgmPCcSncbSJ3r26y1BcaYKmCytTYfwBiTGPTIRET8sWkevPdLp3R35H/BGb9V6W4DCsp32lq7Oxj3ERHx289Kdz+BlFPcjirq+LuwfoK1dkuogwkl7RMRaUQKVvm67q6ErElw/p/VNNElfu0TATYZY4qNMd8aY54HmgFZxpjmIYwtqLRPRKQR8HqdkceUEXCgAK5+Gy7+f0ogLvJ3Oqs9kOH7yMTZZDgdZ6vI98Aia+0vQhGgiAhQq3T3Ihj7NLTQcqzbjptEjDE9rbUbgS98H9XX44EBOEllYIjiExE5vHR37N8h8zqV7oYJf0YinwE9AIwxM3B6Zi3D2RuyBFgSsuhEJLqVFMG/fwcrp0PKELjkOWcRXcKGP+eJ9Kjx6fs4o467gZONMbHAcmC5tfZXoQhQRKLUYaW7D8IZv1Hpbhiq09+Itfalmp8bY7rjJBVNZ4lIcFSWwdw/OsfWdjgRfvEJJKt0N1zVK61ba7cCW4HZwQlHRKJawSp49xbYuQqybobz/6TKqzAXcBIxxsyz1p4VzGBEJEp5vfDdZGcE0rQtXP0OpJ3vdlTih/qMRIYFLQoRiV7Fec6Jg5vnq3Q3BLJzPDw+Zy35RSUktW3GvaPSGZ+ZHLT7a5VKRNyzYjq8/1vwqnQ3FLJzPDwwYwUlFVUAeIpKeGDGCoCgJRJ/d6yHJWNMC2PMv4wxzxtjrnE7HhHxU8lemH4zvHszdEyHX36ltu0h8PictT8mkGolFVU8Pmdt0L5G2CURY8xLxpidxpiVta6PNsasNcasN8bc77t8CTDdWnsLMLbBgxWRuqkohe+ehclD4ftsp3T3pg+19yMErLV4ikqO+Fj+Ua4HIhyns14B/gG8Wn3Btx9lMnAekAcsMsbMAlKAFb6nHZ5uRSR8VJbB0ldh/l+dfR89znQqr5Iy3Y6s0bHW8tX63cccbSS1bRa0r1efJBKScae1dp4xpkety0OA9b72KxhjpgHjcBJKCs7pimE3qhKJepXlkPOakzz2eaD76XDJFEhVYWcoLN26l//7aA3fbdxDcttmTBzSjfdyPJRWeH98TrO4WO4dlR60r1mfJPJl0KI4vmRgW43P84ChwNPAP4wxF3GMvSrGmFuBWwG6d+8ewjBFBICqCsh9A+Y9AcXboNtQGP9PSB2udY8QWLNjH0/M+YFPVxeQ2DKeh8f0ZeLQ7iQ0iWVoaofwrM6y1o4MWhSBx3AQuMmP500BpgBkZWXZUMclErWqKmDZNJj3f1C01Tmqdszf4MSzlTxCYEvhQZ785AdmLcunZUIT7h2Vzo2n96BFwk8/2sdnJgc1adQWjmsiR+IButX4PMV3zW86lEokhKoqYcXb8OVfYO9mSBoEFz0Jvc5V8giBHcWlPP3ZOt5etI0msYbbh5/IbWf1pG3z+AaPxZ9W8P8BPGetLTvGcxKA26y1TwczuBoWAScZY1JxksdVwNV1uYG1djYwOysr65YQxCcSnbxVsOId+PL/YM8G6DoQJr4FaaOUPEJg78FynvlyA//6ZjNVXsvEId256+xedGrd1LWY/BmJdAHWG2M+wFkHWQvsB1oBacAI4AJqVFPVhzFmqu+eicaYPOAha+2Lxpg7gTlALPCStXZVHe+rkYhIsHirYOUMZ+RRuA46D4Cr3oT0C5U8QuBAWSUvfbWJ5+dt5EB5JRMyk7n7nDS6d3D/cFlj7fGXCIwxicCNOMliANAW2IvTBv4D4FVrbWHIogyirKwsu3jxYrfDEIlMXi98/x588RfYvRY69YMR90PviyFGBZLBVlpRxRsLtvLPz9dTeLCcUf0687vz00nr3KpB4zDGLLHWZh3pMb/WRKy1u4EnfB8iEm28Xlg9yxl57PweOvaBy1+BPuOUPEKgssrLu0vz+Nun68gvLuWMXoncMyqdjG5t3Q7tZyJlYb3eNJ0lEgBrYc378MVjULASEtPgspeg7wQljxDwei0frNzOkx//wMbdBxnYrS1PXD6Q03uFb0NKv6azfnyyc676gziL2l2BfGAa8Ii1tjQkEQaZprNE/GAtrP0Qvvhf2LECOvSC4fdD/0sgJtbt6Boday1f/LCLJ+asZVX+PtI6t+Se89M5r29nTBisMdV7OquGZ4B04C5gC3AC8HuczYCT6hOkiIQBa2Hdx/D5/8L2XGiXChOeg/6X6WjaEFm0eQ+Pf7SWhZv30K19M/7flQMZOzCZ2Bj3k4c/6vquGA+caK0t8n3+vTFmAbCeME8ims4SOQZrYf1cZ+ThWQJtT4Bx/4STr1TyCJGVnmKe+HgtX6zdRcdWCfxpfH+uzOpGfJPImias67tjB9AcKKpxrRmwPVgBhYr2iYgcgbWw4TP44lHIWwRtujvnegycCLFxbkfXKG3cdYAnP/mB95dvp02zOO4b3ZsbT+9Bs/jInCasaxJ5DfjIGPN3nP5V3YBfAa8aY86ufpK19rPghSgiQWctbPoSPn8Utn0HrVPg4qcg4xpo0vC7nqNBflEJT89dxztL8khoEsOdI3txy1k9adMsspN1XZPIbb5ff1/r+u2+DwAL6HAAkXC1+StnzWPL19AqCS76q3OiYJMEtyNrlAoPlPHPLzbw2ndbwMJ1p57Ar0b2omOrxvH9rlMSsdamhiqQUNOaiES9Ld84yWPzfGjZBS543DlNMM69lhmN2f7SCp6fv4kX52+kpKKKSwel8OtzTyKlnfu7zIOpriW+Z1lr5x3h+kRr7dSgRhYiKvGVqLN1gbNgvvELaNEJzvwtnHIjxAXvYKJolp3jOazV+t3nnsTeQ+X884sNFB2q4MIBXfjteWn06tSwu8yD6VglvnVNIjuBl4EHrbUVxpi2wHNAprU2LRjBhpqSiESNvMXOyGPDXGjREYbdDVmTIL5x/U/YTdk5Hh6YseJn55gDnJXWkXvPT2dAShsXIguuYO4TycBJIot8i+sP4/TO0hmXIuHA64X1n8B3z8DGz6F5BzjvjzD4FxDfwu3oGp3H56w9YgJJbBnPq5OGuBBRw6vrmki+MWY8sADnkKcXrbW3HftV4UFrItKolRZD7puwcArs2QitusK5D8PgWyChpdvRNTrWWhZu2oOnqOSIjxceKG/giNxTpyRijMkAXsfZXPh74CljzJvAHTU2IIYl7RORRmn3Oidx5L4J5QecY2jPfhD6jNU+jxDYc7Ccd5fkMXXRVjbuOojBKUetLalt9Kw31XU6ay5wn7X2BQBjzOc455yv4PCTB0UkVLxeZ51jwbOw/lOIjYf+l8KQWyF5kNvRNTper+XbjYVMXbiVj1cVUF7l5ZQT2vHE5b3wei0PzVp12JRWs7hY7h2V7mLEDauuSWSwtXZj9Se+M85vNsaMDW5YIvIzpftqTFltcMp0R/6XU2nVspPb0TU6u/aXMX1JHtMWbWVL4SHaNIvjmlO7M3FI98PO84hvEnNYdda9o9JDeqZ5uKnrmshGY8x5wESgo7V2jDEmCzgQkuhEBHav901ZveFMWaUMgZG/d6astLs8qLxey/z1u5m2cCuffF9ApdcyNLU9vzk3jdH9u9A07uetScZnJkdV0qitrmsidwG/Bl4ALvVdLsGZ0jo9uKGJRDGv1+lpteBZp9oqJs6Zshp6KySf4nZ0jU7BvlLeXrSNtxZvI29vCe1bxDPpjFSuHNyNEzuqMOFY6jqddTdwjrV2szHmPt+1NTjt4cOaqrMkIpTth9ypsPA5KFwPLTvDiN87U1atOrsdXaNS5bV8+cNO3lywjc/X7qTKaxnWqwP3je7N+f06k9AkMhsiNrS6JpFWwDbf76uLEuKAsK9nU3WWhLXCDbDwech5Hcr3Q3IWXPIC9B2nKasg8xSV8PaibbyzeBv5xaUktkzg1rN6cmVWN3okai9NXdU1icwD7gceqXHtP4DPgxaRSLTwemHjZ7BginMQVEwT6DcBht4GKUfcHCwBqqzy8tmanUxduJUvf9iFBc48qSP/M6Yv5/TpTFxsZJ3hEU7qmkTuAmYbY24BWhlj1gL7gYuDHplIY1W2H5ZNgwXPQeE6p5/V8Psg6yZo1cXt6BqVbXsO8daibby9eBs795fRqVUCvxrZiyuyutGtvdq/BENdq7O2G2MGA4NxjsbdBiy01npDEZxIo7Jn409TVmX7IGkQXPK8b8qqcbQFDwfllV4+XV3A1IVb+Wr9bgwwIr0TE4d0Z2R6R5po1BFUdT730jodGxf6PkTkWKx1elgteA5+mAMxsb4pq9s1ZRVkm3cfZNqibUxfso3dB8pJatOUX59zEldkdYuqHeQNTYcni4RC2QFYPs1Z79i91umiO/w/4ZSboHVXt6NrNMoqq5izqoBpC7fyzYZCYmMMZ/fuxNVDunNWWkdiY4zbITZ6SiIiwbRnEyx6AZa+BmXF0DUDJjznjD40ZRU063ceYNrCrby7NI+9hypIadeMe85P4/KsbnRurUO2GlLUJBHtE5GQsdY58GnhFFj7oTNl1Xecb8pqMBj9bzgYSiuq+HDldqYu2MbCzXtoEmM4v19nrhrcnTN6JRKjUYcr6nQo1WEvNGaetfasIMcTcjqUSoKm/KBTZbVwCuxaA80TnQqrrEnQOsnt6BqNtTv2M3XhVt7L8VBcUkGPDs25akh3Lh2U0mjOKQ93wTyUqqZh9XitSOTasRJyXoNlU51zPLoOhPHPQL9LdF55kJSUV/H+8nymLtzK0q1FxMfGMKp/FyYO6capqR006ggjUTOdJVIvpftg5XRnrSN/qdN+vffFzsbAbkM1ZRUkq/KLmbZwG9k5HvaXVdKzYwsevKgPlwxKoX0L7dwPR0oiIkdjLWz91kkcq96DyhLo1BdGPQonXwktOrgdYUTKzvEc1jr9rrN7YYFpC7eyLK+Y+CYxXDSgKxOHdGdwj3YYJeiwpiQiUtv+Alj2prMpsHA9xLeCgVdC5vXOoU/6oRaw7BwPD8xY8eMhTp6iEu6fsQKA9M6teGhMXyZkJtO2uUYdkUJJRASgqtJpub70NfjhI7BV0P00OOO30G88xKsxXzA89uGaw04BrJbYMoGP7j5To44IVJ8kor9tiXyFG5wRR+6bcGCHsynwtF9B5nXQMc3t6BqFg2WVzFm1g/dyPOzYV3rE5xQeKFMCiVD1SSJfBi0KkYZUfghWz3JGHVu+AhMDJ53vJI60URAb53aEEa+yysvXGwp5b2kec1YVUFJRRUq7ZrRKaML+ssqfPV9tSSJXwEnEWjsymIGIhJS1sD3XSRwrpju7ydulwtn/DRlXa19HEFhrWZW/jxlLPcxals/uA2W0aRbHhEHJTMhMJuuEdszMzT9sTQSgWVws944K+3Pt5Ci0JiKNW8leWP4OLH0VClZAk6bO2eSDrocThkGMOrrWV97eQ8zMzee9HA/rdx4gLtbpXzUhM4WRvTsedkJg9VnkNauz7h2VHtVnlEe6gHeshwNjTE/gv4A21trL/HmNdqxHAa8XNs9zRh2rZ0NVmbMhMPM6GHA5NGvrdoQRr7ikgg9XbOe9HA8LNu0BYHCPdozPTOaiAV1VXdXIhGTHunFWwW6w1r4S4OtfwjnMaqe1tn+N66OBvwGxwAvW2seOdg9r7UbgZmPM9EBikEam2OMskOe8BkVboGkbZ8Qx6DoniUi9lFd6+fKHXbyXk8enq3dSXumlZ2ILfndeGuMzk3XIU5Sqz5qI9TU1fCXAW7wC/AN4tfqCMSYWmAycB+QBi4wxs3ASyqO1Xj/JWrszwK8tjUVluVOSu/RV2DAXrBdSz3LWOvpcDHFasK0Pay1LtxaRnePh/eX57D1UQYcW8Vw9pDsTMpM5OaWNqqqiXH3XROKMMV8AiwEvgLX2P/15obV2njGmR63LQ4D1vhEGxphpwDhr7aPU4wheY8ytwK0A3bt3D/Q2Ek52rXUSx7JpcGg3tEpy9nRkXgvtU92OLuJt2n2Q7BwP2bkethQeIqFJDOf368IlmcmccVKiziSXH9U3ify11uf1XWBJxjlyt1oeMPRoTzbGdAAeATKNMQ/4ks3PWGunAFPAWROpZ4zilrIDTvuRpa9C3kKIaQJpo2HQDdDrHKcFuwRsz8Fy3l/uLJDnbC3CGDitZwfuHNmL0f270KqpSp/l5wJKIsaY26y1z+GMDmr/UJ5X76j8ZK0tBG7357k6TyRCWQt5i2Hpv5wEUn4AEtPgvD/BwKugZSe3I4xopRVVzF29k/dy8vhi7S4qvZbeXVrxwAW9GZuRRNc2mg6UYwt0JPKd79f3gxWIjwfoVuPzFN+1erPWzgZmZ2Vl3RKM+0mIHdztTFXlvOac1RHX3Gm1Puh66DZE/avqweu1LNi0h+wcDx+s2M7+sko6t05g0hmpTMhMpk/X1m6HKBEkoCRirV3m++1KYI9vkd0A7esZzyLgJGNMKk7yuAq4up73lEjhrYINnzujjrUfgrfCORlwzNPQ/xJIaOV2hBGhdpfc6n0YPxTs570cDzNzPOQXl9IiPpbR/bsyITOZ007soPPIJSD12idijJlrrT3naJ8f57VTgRFAIlAAPGStfdEYcyHwFE5F1kvW2kcCDvDwr1c9nXXLunXrgnFLCZa9W37qX7UvD5p3gJOvckpzO/VxO7qIUrtLLkBcrKFTqwQ8RaXExhjOPCmRCZnJnNe3M83jtd9Yju9Y+0Tqm0QOOyLXGDPfWntmwDdsANpsGCYqSmHN+8501UZfG7YTz3amq9IvhCbarBaIYY99hqeo5GfX42IND1zQhzEDk3SkrNRZqI7HBVhujPkbTjPG4cDyet4vZLSwHiZ2rHSqq5a/BaVF0KY7jHjA6V/VtttxXy5HVlHlZd4Pu46YQAAqqyyTzlDpswSfX0nEGHMvsBTIsdbuqb5urb3T98O5D/Cpb/E6LGlh3UUlRbDyXWfUkZ/z09Gyg66D1BHqXxUgay1LtuwlO9fDv5dvZ++hCoxxCtpqU5dcCRV/RyKjgPuAdsaYPJyEsgiYWf3DOUTxSaSqqoD1c2HZVGeRvKoMOvWD0X+Bk6+A5vWtwYhe6wr2k53rYWZuPnl7S2gaF8O5fTozPiOZokPl/PfMVeqSKw3GryRirT0XwBhzApAJDALOAv7b15bkJmvtoZBFKZHBWtix3CnNXfEOHNzlLJJn3eScSZ6UqdLcAG0vLmFWbj7Zufms3r6PGAPDeiXym3PTGNW/Cy0Tfvqn3CQ2Rl1ypcHUd2E9EXgTWGyt/X3QogoBVWeF0L7tTtJYNhV2fu9MV6WNhoET4aTzdMhTgKo75WbnOp1yrYWB3doyPiOJi0/WArk0nJBVZ/lungZ8YK2NiBVrVWcFSfkhWPuBU5a78XOn8WHKYGcXeb9LNF0VoNKKKj5fs5PsXA+fr9lFeZWX1MQWjMtIYlxGMqmJOutdGl4oq7MAtgJdg3AfCXdeL2z9xhlxrJoJ5fuhTTen8eHAiZAYEf+PCDtVXst3GwvJzvHw0cod7C+rpGOrBK499QTGZyYxIFmdciV8+VudVQTk+D6W+n5dba31AtcAG0IVYLCoxLceCjc46xzLp0HRVohvCX3HO6MOnQ4YkOqjZLNzPMxenk/BvjJaJjRhVL8ujM9M4rSeHWiiTrkSAfyazjLGDMNZUM8EMoB+OI0XS4AE4Apr7b9DF2bwaDrLTyV7nYaHy6bBtgWAgRNHOiOO3hdBvKZVArG18BAzc50W6xt2HSQu1jA8rRPjM5M4t09nmsapE7GEn3pPZ1lrvwa+rnHDJkBfoBOwwlpbEIxAxWU/luW+6SvLLYeOveHcPzhlua2T3I4wIhUeKOPfK7aTneNh6dYiAIaktmfSGak6SlYiXqANGCsJ493pUgfWwvZlP5XlHtoNzRMh62ZnuqrrQJXlBuBQeSUfryogO9fD/HW7qfK1WL9vtNNiPVmb/6SROG4SMcZcX8d75lprlWDC3b7tsOJtJ3lUl+WmX+BMV/U6V2W5x3C0LrkVVV6+Wreb7FwPH68qoKSiiqQ2TbnlzJ6Mz0yidxe1WJfGx5+RyEj8O7HQ+J5XRBiOUrSwjlOWu+bfTnXVj2W5Q+CiJ6HfBJXl+qF2l1xPUQn3vbuc6Uu28f32/ew5WE6bZnGMz0xmfEYSg3u0J0Yt1qUR8yeJbMZJDv78S6hOImEnantnVZfl5k6F76vLcrvDmb9zRh0dTnQ7wojy+Jy1h7UUASir9PLV+kIuOrkr4zOSGZ7WkfgmqqyS6ODvmoi//5XSf7nCReEGZ8Sx7C0orlGWmzERup+ustwA7CguPWqXXANMvnpQwwYkEgaOm0SstX9oiEAkCA7t+aksN28hmBjoOQLO+R9fWW5ztyOMOMWHKvhw5U+tR45GXXIlWulYs0hXWQ7rP3ESxw8f+cpy+8B5f4QBl6ssNwClFVXMXb2TmbkevljrtB7pmdiCu89Jo3l8LE9+8oO65Ir4KIlEImth20LnYKdVM5yNgSrLrZfKKi/fbChkZm4+c1bt4EBZJZ1aJXDdaScwPiOZ/smtf2w90rFVgrrkivhETRJpFNVZu9c7ZbnL34K9m6FJM2ea6uQrnd3kKsutE2sty/KKyc7x8P7y7ew+UEarhCZcOKAL4zKSObVnB2KPUFk1PjNZSUPEp95dfCNNxLU9ObgbVs5w+lZ5lgAGeg6Hk6+CPhdDQiu3I4w4G3YdYGZuPrNyPWwuPER8bAxn93Zaj4xI76TWIyK1hLqLrwRbdZv15W/D+k/BVkHnAXD+n6H/pVrnCEDBvlJmL8tnZm4+KzzFGAOnn9iBO0b0YlT/LrRpplGcSCCURMKFtwo2z3cSx/eznP0crZLg9Dud6arO/dyOMOIUl1QwZ+UOsnM9fLuxEGvh5JQ2PHhRH8YMTKJz66ZuhygS8ZRE3LZjpbPGsWI67M+H+FbQdxwMvNLXZl1TK3VRfajTzNx8Plu7k/JKLz06NOc/zj6JsRlJnNixpdshijQqSiJuKPbAyunOqKNgJcQ0gV7nwahHnP5VcdpzUBdVXsu3GwqZmfvToU6JLRO4Zmh3xmckc3KKDnUSCRUlkYZSug9Wz3YWyDfNB6xznOyFTzh9q1okuh1hRLHWssJTTHZOPrOX57Nrv3Oo0+j+XRiXoUOdRBqKkkgoVVXAhs+cjYBrP4DKUmiXCsPvc87nUN+qOtu0+yAzcz3Mys1n4+6DxMfGMLJ3R8ZlJHN2b1VWiTS0qEkiDbZPxFrwLHXWOVa+65zP0awdZF7rLJCnDNZGwDraub+U2cu2MyvXw7I8p7Lq1NQO3Da8J6P7daVNc1VWibhF+0SCZc8m51Cn5W9B4XqITXDWN06+0jmfo4lOr6uLfaVOZdXM3Hy+2bAbr4V+Sa0Zn5HMxQO70rWN1o1EGor2iYRKdcPD5W/5ziEHepwJw34NfcZCs7auhhdpyiqr+HzNLmbmepi7xqms6t6+OXeO7MXYjCR6ddLGSpFwoyRSVxWlTqPD5W/Duo/BW+GcQ37OQ07Dw7bd3I4wolR5LQs2FTIzJ58PVm5nf2kliS3juXpId8ZlJJHRra0qq0TCmJKIP6yFLV/7Gh7OhLJiaNkZht7mTFd1GaB1jmOofZzsPeencVLnVmTneJi9PJ+CfWW0iI9lVH+nZ9WwE1VZJRIplET8NfNOOLAT+oxxKqt6jtBGQD8c6TjZ3769DAvExRqGp3Xivy9O4pzenWkWr++nSKRREvGHMXDVG9CuB8S3cDuaiPLYh2t+dpysBdo2i+OLe0fQtrkKDkQimZKIv9S7ym/7SyuYs6qAmbkeduwrPeJziksqlEBEGgElEQmKssoqvly7i5m5+Xy6uoCySi/d2jejVUIT9pdV/uz5Ok5WpHFQEpGAeb2WBZv2MDPXwwcrtrOvtJIOLeK5anA3xmYkM6h7W2bm5h+2JgI6TlakMYnoJGKMGQ9cBLQGXrTWfuxuRI2ftZZV+fuYmeth9rLt7NhX6lRW9evC2IwkhvVKJK5GZVX1CYA6TlakcXJtx7ox5iXgYmCntbZ/jeujgb8BscAL1trH/LhXO+AJa+3Nx3tuxJ1sGCa2FB5kZm4+M3M9bNh1kCYxhhHpTs+qc/uoskqkMQvXHeuvAP8AXq2+YIyJBSYD5wF5wCJjzCychPJorddPstbu9P3+Qd/rJIh27S/j38vzyc7NJ3dbEQBDUttz8xk9uaB/F9q10MK4SLRzLYlYa+cZY3rUujwEWG+t3QhgjJkGjLPWPoozajmMcbYyPwZ8aK1dGuKQo8L+0go+XlVAdq6Hr9c7Pav6dG3NAxf0ZszAJC2Ii8hhwm1NJBnYVuPzPGDoMZ5/F3Au0MYY08ta++yRnmSMuRW4FaB79+5BCrXx+LGyalk+n37/U2XVHSOcnlVpndWzSkSOLNySSJ1Ya58GnvbjeVOAKeCsiYQ6rkhQXVk1a5mHD1bsoLikgvYt4rlycDfG+Sqr1LNKRI4n3JKIB6jZwTDFd63eGuw8kTBWXVk1a1k+s3Lz2bGvlOa+yqpxR6isEhE5nnBLIouAk4wxqTjJ4yrg6mDc2Fo7G5idlZV1SzDuF0m2FB5kVm4+M5fls37ngR8rq35/UR/OU2WViNSDa0nEGDMVGAEkGmPygIestS8aY+4E5uBUZL1krV0VpK/XaEcitbvk3jsqnTNOSuT9ZU7iyNlaBDiVVY9M6M+F/buqskpEgkInG0a42l1yAWKM073e4lRWjctIYszAJJJVWSUiAQjXfSISBP/30c+75HottExowow7TldllYiEVNQkkcY0nVWzsiq/+Mhdcg+WVSqBiEjIRU0SifSFdWst32/fx8zcwyurmsXF/mwkAuqSKyINI2qSSKTaWniImbmewyqrhqc5lVXn9unEx6sK1CVXRFwTNUkkkqazdh8o49/Lt5Od6/mpsqpHe/48vj8XDTi8skpdckXETarOChMHyir5eNUOsnPz+Xr9bqq8lt5dWjE+M1mVVSLiKlVnhanySi9f/rCLmbkePl1dQGmFl+S2zbjtrJ6My0gmvYsWxkUkvCmJNDCv17Jw8x5m5ubzwYrtP/asuvyUbozPTGJQ93bqWSUiESNqkoibayLVlVWzcvOZtSyf7cVOZdX5fTszLjOZM9SzSkQilNZEQmjbHl9lVW4+62pUVo3NSOK8vp1pHh81OVxEIpjWRBpQdWXVzFwPS49TWSUiEumURIKgurJqZm4+X9WorLpvdG/GDOxKSrvmbocoIhISSiJ+OFKX3AsHdGXeD7vIVmWViESxqFkTqbGwfsu6dev8ft2RuuTGxhjiYw0lFV7aNY/j4pOTGJfhVFbFxKiySkQaF62JEHjvrMfnrP1Zb6oqr4XYGF6+cTBnnKTKKhGJXlGTRAKVX1RyxOulFVWM7N2pgaMREQkv+i/0cRytG6665IqIKIkc172j0mkWd/gZ5OqSKyLiiJrprEB3rKtLrojI0UVNdVa1cO3iKyISro5VnaXpLBERCZiSiIiIBExJREREAqYkIiIiAVMSERGRgEVddZYxZhewBWgDFNfx5f6+xp/nHe85x3r8aI8lAruPG517AvmeN9S9Q/l+8Pe5gfydH+sxvR8a9h6N+f1wgrW24xEfsdZG5QcwJVSv8ed5x3vOsR4/2mPAYre/r8H+njfUvUP5fgjGe0Lvh4a/d13vEa3vh2iezpodwtf487zjPedYjwcSezgIZdz1vXco3w/+PjfQv3O9H0Jz77reIyrfD1E3ndWYGWMW26NsCJLoo/eD1BSq90M0j0QaoyluByBhRe8HqSkk7weNREREJGAaiYiISMCUREREJGBKIiIiEjAlkShhjBlvjHneGPOWMeZ8t+MRdxljehpjXjTGTHc7FnGHMaaFMeZfvp8L1wR6HyWRCGCMeckYs9MYs7LW9dHGmLXGmPXGmPuPdQ9rbba19hbgduDKUMYroRWk98NGa+3NoY1UGlod3xuXANN9PxfGBvo1lUQiwyvA6JoXjDGxwGTgAqAvMNEY09cYM8AY836tj041Xvqg73USuV4heO8HaVxewc/3BpACbPM9rSrQLxg1x+NGMmvtPGNMj1qXhwDrrbUbAYwx04Bx1tpHgYtr38MYY4DHgA+ttUtDHLKEUDDeD9I41eW9AeThJJJc6jGg0EgkciXz0/8iwHlDHOvg97uAc4HLjDG3hzIwcUWd3g/GmA7GmGeBTGPMA6EOTlx1tPfGDOBSY8wz1KNVikYiUcJa+zTwtNtxSHiw1hbirI9JlLLWHgRuqu99NBKJXB6gW43PU3zXJDrp/SBHE9L3hpJI5FoEnGSMSTXGxANXAbNcjknco/eDHE1I3xtKIhHAGDMV+BZIN8bkGWNuttZWAncCc4DVwNvW2lVuxikNQ+8HORo33htqwCgiIgHTSERERAKmJCIiIgFTEhERkYApiYiISMCUREREJGBKIiIiEjAlERERCZiSiIiIBExJREREAqYkIuIiY8xqY8wBY0y57+OA76OP27GJ+ENtT0TCgDHmRWCjtfYRt2MRqQuNRETCw8nAyuM+SyTMKImIuMwYE4Nz9rWSiEQcJRER93XH+be40e1AROpKSUTEfa2Bg0C824GI1JWSiIj7VgPLgL3GmN5uByNSF6rOEhGRgGkkIiIiAVMSERGRgCmJiIhIwJREREQkYEoiIiISMCUREREJmJKIiIgETElEREQCpiQiIiIB+/99V0RrTBKXHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 计算两个酉矩阵之间的 L1 谱距离，即 (5) 中定义的误差\n",
    "def calculate_error(U1, U2):\n",
    "    return np.abs(linalg.eig(U1 - U2)[0]).max()\n",
    "\n",
    "# 封装计算误差的函数，自由参数为每个时间块的长度 tau 和 product formula 的阶数 order\n",
    "def calculate_total_error(tau, order=1):\n",
    "    # 注意因为在电路中的多 Pauli 旋转门时会引入一个额外的全局相位，需要在计算误差时扣除\n",
    "    # 补充：该全局相位不会对实际的量子态产生任何可观测的影响，只是计算理论误差时需要扣除\n",
    "    h_2 = Hamiltonian([[1, 'Z0, Z1'], [1, 'X0'], [1, 'X1']])\n",
    "    cir = UAnsatz(2)\n",
    "    # 应为总时长为 1，故 steps = int(1/tau)，注意传入的 tau 需要可以整除 1\n",
    "    construct_trotter_circuit(cir, h_2, tau=tau, steps=int(1/tau), order=order)\n",
    "    cir_U = cir.U.numpy()\n",
    "    U_evolve =  np.exp(1j) * linalg.expm(-1 * 1j * h_2.construct_h_matrix()) #理论上的时间演化算符加上一个全局相位\n",
    "    return calculate_error(cir_U, U_evolve)\n",
    "\n",
    "# 不同的参数 tau，注意它们需要可以整除 1\n",
    "taus = np.array([0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1])\n",
    "errors = []\n",
    "\n",
    "# 记录对应不同 tau 的总误差\n",
    "for tau in taus:\n",
    "    errors.append(calculate_total_error(tau))    \n",
    "\n",
    "# 可视化结果\n",
    "plt.loglog(taus, errors, 'o-', label='error')\n",
    "plt.loglog(taus, (3 * taus**2 * (1/taus) * np.exp(3 * taus)), label='bound') # 按照 (10) 计算的一阶误差上界\n",
    "plt.legend()\n",
    "plt.xlabel(r'$\\tau$', fontsize=12)\n",
    "plt.ylabel(r'$\\Vert U_{\\rm cir} - \\exp(-iHt) \\Vert$', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fdbf168",
   "metadata": {},
   "source": [
    "下面，我们将 `order` 设置为 2，并计算二阶 product formula 的误差："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a3db610c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAENCAYAAADOhVhvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAw4ElEQVR4nO3deXxU9b3/8dc3C0nY17AlYREIO4QEcENAVEBkE6z7UrVurbe21dtqbW3vrdVbrb/Wq7fVum+gBkVcWlRccKECWdhBkXUSICSQEELWme/vjxMghEC2mZyZ5P18PPKAOXPmzCfhkM98t8/XWGsRERFpiDC3AxARkdClJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg0W4XYATa1r1662b9++bochIhJS0tLScq213aofb3FJpG/fvqxevdrtMEREQooxZmdNx9WdJSIiDaYkIiIiDaYkIiIiDdbixkRqUl5ejsfjoaSkxO1QXBUdHU1cXByRkZFuhyIiIUJJBPB4PLRr146+fftijHE7HFdYa8nLy8Pj8dCvXz+3wxGREKHuLKCkpIQuXbq02AQCYIyhS5cuLb41JiL1oyRSqSUnkKP0MxBppnxeyFwAPp/fL63uLBGR5qwoFxbdBNs+g5hOkDjNr5dXEmmAxRlZPLJ0C9n5xfTqGMM9UxOZk9Q7oO/p9XoJDw8/5eOaWGux1hIWpganSIvkWQ1vXA9F+2HWE35PIKDurHpbnJHFvW+tIyu/GAtk5Rdz71vrWJyR1ajrvvLKK4wbN47Ro0dz66234vV6adu2Lb/4xS8YNWoUK1asOOnxY489xvDhwxk+fDh/+ctfANixYweJiYlcd911DB8+nN27dzf+mxaR0GItrPwHPDcNwsLgpg9hzLUBeSu1RKr5/bsb2Jh96JTPZ+zKp8x7Yr9icbmX/0xdy4KVu2p8zdBe7Xlg5rBTXnPTpk28/vrrfPXVV0RGRnLHHXfw6quvUlRUxPjx4/nzn/8McMLjtLQ0nn/+eb755hustYwfP56JEyfSqVMnvvvuO1588UXOPPPMBvwERCSklR2B9+6Cta/DwItg7lPQunPA3k5JpJ6qJ5DajtfFsmXLSEtLY+zYsQAUFxcTGxtLeHg48+bNO3Ze1cdffvklc+fOpU2bNgBceumlfPHFF8yaNYs+ffoogYi0RHnfw+vXQs5GmPxrmHC30xIJICWRak7XYgA45+FPyMovPul4744xvH7rWQ16T2st119/PQ899NAJxx999NETxj2io6NrHQcBjiUWEWlBNr8Pb98GYeFwTSoMuKBJ3lZjIvV0z9REYiJP/EUeExnOPVMTG3zNKVOmkJqaSk5ODgAHDhxg584aC2YeM2HCBBYvXsyRI0coKiri7bffZsKECQ2OQURClLcCPnoAFl4FXc6AW5c3WQIBtUTq7egsLH/Ozho6dCh/+MMfuOiii/D5fERGRvLkk0+e9jVjxozhhhtuYNy4cQDcfPPNJCUlsWPHjgbHISIh5vB+SP0h7PgCkn8I0/8HIqKaNARjrW3SN3RbSkqKrb6fyKZNmxgyZIhLEQUX/SxEQsTulc703eIDcMn/g9FXBfTtjDFp1tqU6sfVEhERCSXWwsqnYel90CEObv4YeoxwLRwlERGRUFFWBEv+A9anwqDpMPfvENPR1ZCUREREQkHuVnj9GsjdAuf/Bs79ecCn79aFkoiISLDbuAQW3wERreCat+CMyW5HdIySiIhIsPJWwLLfwdf/C71T4AcvOuMgQURJREQkGBXuc6bv7vwKxt4MU//Y5NN368L9DjUBnMKJw4cPb7L3+93vfsejjz7aZO8nIvWwcwU8dR5kpcPcp2HGn4MygYBaIiIiwcNa+Pff4KPfQMcEuPYt6H76UkxuU0skiFRUVHD11VczZMgQ5s+fz5EjR1i2bBlJSUmMGDGCG2+8kdLSUgD69u1Lbm4uAKtXr2bSpEmA08K48cYbmTRpEv379+fxxx8/dv0HH3yQQYMGce6557Jly5Ym//5E5DRKC53uq6X3wqBpcMtnQZ9AQC2Rk/3zV7B3nX+v2WMETH+41tO2bNnCs88+yznnnMONN97IY489xlNPPcWyZcsYNGgQ1113HX/729+46667TnudzZs38+mnn1JYWEhiYiK33347a9euZeHChWRmZlJRUcGYMWNITk720zcoIo2yf4tTfTfvO7jg93DOTyFEtqtWSySIxMfHc8455wBwzTXXsGzZMvr168egQYMAuP7661m+fHmt15kxYwZRUVF07dqV2NhY9u3bxxdffMHcuXNp3bo17du3Z9asWQH9XkSkjja8Df843ylfct07cO5dIZNAQC2Rk9WhxRAoptqN07FjR/Ly8mo8NyIiAp/P2cOkpKTkhOeioo4PwIWHh1NRUeHnSEWk0bzlTvXdfz8JceOc6bvte7kdVb2pJRJEdu3axYoVKwB47bXXSElJYceOHWzduhWAl19+mYkTJwLOmEhaWhoAixYtqvXa5513HosXL6a4uJjCwkLefffdAH0XIlKrwr3w4kwngYy7FW54PyQTCCiJBJXExESefPJJhgwZwsGDB/nZz37G888/z2WXXcaIESMICwvjtttuA+CBBx7gpz/9KSkpKXXaqGrMmDFcfvnljBo1iunTpx/bRVFEmtiOr+DvE2DPGpj3LFz8J2cleohSKXhU/rwq/SxEAsRaWPGE04XVuR9c/grEhs7/NZWCFxFxS8kheOfHsGkJDJkJs/8Potu7HZVfKImIiARSziZn+u6BbXDRH+Csn4TU7KvahHQSMcbMAWYA7YFnrbUfuhuRiEgV61JhyZ3Qqi1cvwT6nut2RH4XdAPrxpjnjDE5xpj11Y5PM8ZsMcZsNcb8CsBau9ha+yPgNuDyxrxvSxsbqol+BiJ+UlEG//wlLLoJeoyEW5c3ywQCQZhEgBeAaVUPGGPCgSeB6cBQ4EpjzNAqp9xf+XyDREdHk5eX16J/iVprycvLIzo62u1QRELboWx48RL45u9w5h1ww3vQvqfbUQVM0HVnWWuXG2P6Vjs8Dthqrd0GYIxZCMw2xmwCHgb+aa1NP9U1jTG3ALcAJCQknPR8XFwcHo+H/fv3++ebCFHR0dHExQXXXgUiIWX7cki9EcqOwPznYfilbkcUcEGXRE6hN7C7ymMPMB64E7gA6GCMGWCt/XtNL7bWPg08Dc4U3+rPR0ZG0q9fP78HLSIthLcCvn4cPvlv6DLAWTzYLdHtqJpEqCSRGllrHwcer/VEEZFAsBa+XeqUbs/9FobOgdlPQFQ7tyNrMqGSRLKA+CqP4yqPiYi4Y88a+PB+pwurywC44jVIvLhZTd+ti1BJIquAgcaYfjjJ4wrgKndDEpEWqSDL6bZasxBiOsH0RyDlhxAe6XZkrgi6JGKMWQBMAroaYzzAA9baZ40xPwGWAuHAc9baDS6GKSItTWkhfPkXp3SJtXDOf8CEX0B0B7cjc1XQJRFr7ZWnOP4B8EEThyMiLZ23AjJegk//CEX7YcRlcP5voFMftyMLCkGXREREgoK18N2H8OFvIHcLJJwNV70OvbUjaFVKIiIi1e1ZWzlo/jl0PgMufxUGz2hxg+Z1oSQiInJUQRZ88gdYs6By0PxPkHJjix00rwslERGR0kL46q/w9RNgvc6g+bk/h5iObkcW9JRERKTl8lZAxsuVg+Y5MHw+TPmtBs3rQUlERFoea+G7j5yV5vs3Q8JZcOUCiDtp4z6phZKIiLQse9Y6yWPbZ9C5v7NN7eBLNGjeQEoiItIyHMp2Bs0zX3PGOqb9jzNoHtHK7chCmpKIiDRvpYXw1ePw9f86g+Zn3+msNNeguV8oiYhI8+StgMxX4JMHKwfN51UOmvd1O7JmRUlERJoXa2Hrx85K8/2bIP5MDZoHkJKIiDQfe9c5yWPbp86g+Q9ehiEzNWgeQEoiIhL6Du2pHDR/tXLQ/GFIuUmD5k1ASUREQlfpYWdb2q//F3wVcNaP4by7nZIl0iRqTSLGmOvqec1Ma+3aBsYjIlI7nxcyXoFPH4TD+2DYpc6geed+bkfW4tSlJTIZsHU4z1Selw8oiYhIYHz3sbNYMGcjxI93KuzGj3U7qharLklkB05yqMvI1NEkIiLiX3vXO8nj+0+gUz/4wUswZJYGzV1W1zGRuv4r6V9TRPzr0B749A+Q8aqzFe3Uh2DszRo0DxK1JhFr7e+bIhARkROUHnYGzL9+HLzlGjQPUnWenWWMiQWmAqOAjjjdVmuAj6y1ewMRnIi0QD6vM1X3kwfh8F4YNhemPKBB8yAVVtsJxpghxphUYBNwLRAJ7K3881pggzEm1RgzNKCRikjzt/Vj+PsEWHIndEyAmz6Cy15QAglidWmJvAA8AlxtrS2t/qQxJgqYBTwLnOXX6ESkZdi3wVlp/v0yp7bVZS/C0NkaNA8BdRkTGX/078aYMGutr9rzpcCblV8iInVXuPf4SvOo9jD1j5WD5lFuRyZ1VN8V6/uAboEIRERakLIiZ9D8q786g+bjb3cGzVt3djsyqac6JRFjzChgIxB9iud3WWsT/BmYiDRDh/fDyqdh1TNQfACGzoELHnCKJUpIqmtL5D0gFggzxiwAMnFmZmXiDM53CERwItJM5H4HK56AzAXgLYPEi+Hcn2mleTNQpyRirY03xnQFdgFf4EzzvRQYjtM6+XvAIhSR0GQt7FrhdFtt+QDCo2D0Vc56j64D3Y5O/KTOYyLW2lxjzAhr7fdHjxljDBBjrT0SkOhEJPT4vLBpiZM8stIgpjNM/CWM/RG01ZBqc1PXMZGncLquMowxe44mDWutBZRARMQZLM941em2yt/pjHPM+DOMugpatXY7OgmQurZESoDLgN8BnY0xW6lMKpV/ZlprcwIQn4gEu8J9xwfLS/IhbhxMfdAZ9wgLdzs6CbC6jon89OjfjTHdcVavbwNSgNuBBEB3i0hLkrPZaXWsfd2ZpjvkEjjrTkgYX/trpdmo986G1tp9xhgv8KS1NhugctBdRJo7a2HHl854x3dLISIakq51Bsu7nOF2dOICv2yPa63N9cd1RCRIeStg0ztO8sjOgNZdYdJ9MPYmaKPPkC1ZXQfW+1hrdwY6GBEJMqWFzja0K/4PCnZBlwFwyV9g1BUQGeN2dBIE6toS2W6MKcRZtb4eiAFSjDEfa3qvSDN0aA+sfApWPwclBZBwFkx/GAZNh7Bai39LC1LXJNIZGF35lQR8D6TiLBXZCKyy1t4ciABFpAnlbHK6rNa+AdYLQ2Y6g+VaWS6nUGsSMcb0t9ZuAz6r/Dp6vBUwAiepjApQfCISaNbC9uVO8tj6EUTEQPINcNYdqmkltapLS+QToC+AMeYtnJpZa3DWhqQBaQGLTkQCx1sOGxY728/uXQttusHk+53BclXTlTqqy34ifas8fA+n1XEXMNIYEw6sBdZaa38ciABFxM9KDkH6S/Dvv8EhD3QdBDMfh5GXQ2SNhbpFTqleU3yttc9VfWyMScBJKq50Zxlj2gD/B5QBn1lrX3UjDpGQcCjbSRxpL0DpIehzjlOWZOBFGiyXBmvUOhFr7S6cyr7v+iccMMY8B1wC5Fhrh1c5Pg34K87K+GestQ/jVBJOtda+a4x5HVASEalu73pnZfm6N8H6nG1nz7oT4pLdjkyagQYnEWPMcmvtef4MptILwBPAS1XeKxx4ErgQ8ACrjDFLgDhgXeVp3gDEIhKarIVtnzqD5d9/ApFtnG1nz7zd2cNcWozFGVk8snQL2fnF9OoYwz1TE5mT1Ntv129MS+Qcv0VRhbV2uTGmb7XD44CtlbPEMMYsBGbjJJQ4jm+OVSNjzC3ALQAJCdqAUZqxkkNOLatVz8D+zdC2O0z5LST/UIPlLdDijCzufWsdxeXOZ+ys/GLufcv53O2vROKXsidNoDewu8pjDzAeeBx4whgzg9N0qVlrnwaeBkhJSbEBjFPEHXvXwapnnfUd5UXQczTMfhJGXAYRUW5HJy55ZOmWYwnkqOJyL48s3dLikkiNrLVFwA/djkPEFRWlsHGJ0+rY/W+nGOLwec4U3d4a7xDIzi+u1/GGCJUkkgXEV3kcV3lMpOU5uBPSnof0l+FIrrMg8KIHna1n1WUlQMaugzzzxXZO1e3Sq6P/6p6FShJZBQw0xvTDSR5XAFe5G5JIE/L54PtlTqvj26VgjFPHauxN0H+ypugKXp/lo417eeaL7azeeZB20RGcP7gbX2/No6TCd+y8mMhw7pma6Lf3bUwSMX6LoupFjVkATAK6GmM8wAPW2meNMT8BluJM8X3OWrshEO8vElSK8iDjZacQYv5OaBML593tlCXpEOd2dBIEikorSE3z8NxX29mZd4T4zjE8MHMoP0iJp01URMBnZxlnm/QGvNCYT621k/0WSRNJSUmxq1evdjsMkVOzFjyrnVbHhrfBWwp9zoWxN8LgmRDRyu0IJQjsO1TCC1/v4LVvdlFQXM6YhI78aEJ/LhrWg/Aw/3/GN8akWWtTqh9vcEskFBOISFArK3IWBK56xplt1aodjLnO6bKKHeJ2dBIkNmYf4pkvtvHu2my8Psu04T246dz+JPfp5Eo8oTImItJ87f8WVj8LmQugtABih8GMx2DkDyCqndvRSRDw+Syff7ufZ77cxldb82jdKpyrx/fhxnP6kdCltaux1aUU/H8AT1lrS09zThRwq7X2cX8GJ9Jsecth8/tOq2PHFxAWCcPmOKvK48c7A+fS4pWUe1mckcUzX25na85herSP5lfTB3PluAQ6xES6HR5Qt5ZID2CrMeYD4HNgC1AItAMG4QyCT6dKmRIROYVD2ZD2olME8fBe6JDgrChPug7adnM7OgkSeYdLefnfO3l5xU7yisoY2rM9/+/yUcwY0YtWEcE1E68upeDvM8Y8BtwA3ISzEVVH4CBOGfgPgPustXmBC1MkhFkL2z93Wh2bP3CKIA64AMb+FQZeCGHhbkcoQWJrzmGe/XIbi9KzKKvwcf7gWG6e0I+z+nfBBGnrtE5jItbaXODRyi8RqYvifFizwClHkvcdxHSGs3/i1LHq3M/t6CRIWGtZ8X0ez3y5nU825xAVEca8MXHcdG4/BsS2dTu8WmlgXcTfsjOdgfK1b0JFMcSNhblPwdA52vRJjimr8PH+umz+sXw7G/ccokubVvzsgkFcc2YCXdqGTr2zeiWRyn3V78dZLd4TyAYWAg9aa0v8H55IiCgvcdZ0rHoGslZDZGtndtXYm6CnK3u2SZAqOFLOayt38eLXO9h7qIQBsW15+NIRzEnqTXRk6HVt1rcl8jcgEbgT2An0Ae7DqbJ7o39DEwkBB7Y5q8kzXoXiA9BlIEz7Hxh1BcR0dDs6CSK78o7w3FfbeWP1bo6UeTlnQBcemjeCiQO7ERaAxYFNpb5JZA5whrU2v/LxRmPMN8BWlESkpfCWw7f/cmZYbf0YTDgMucSZntt3gqbnygnSdh7gmS+2s3TDXsLDDDNH9eLmc/sztFd7t0Pzi/omkb1AayC/yrEYYI+/AhIJWnnfQ/pLkPkaFOVAu54w6V4Ycz207+l2dOKi6vWpfnHhQKJbRfCPL7aRsSufDjGR3DbxDK4/uy/d2zevcbH6JpGXgX8ZY/4XZ2OoeODHwEvGmPOPnmSt/cR/IYq4qLzY2bMj/SXY+aXT6hg0zSlHMuACCNfclJaupt0Df/HmWiyQ0Lk1v581jPnJcbSJap73Sn2/q1sr/7yv2vHbKr8ALNC/MUGJuG7vOidxrH0dSgqgUz+Y8oCzZ0e7Hm5HJ0Gkpt0DLdC5TSs+vXtSQIohBpN6JRFrrSa3S/NVcgjWpzrJIzsDwqNg6Gyn1dHnHO3ZIScoKfeydMNesk6xS+DBorJmn0Cg/lN8z7PWLq/h+JXW2gX+C0ukiVgLu79xEseGt6H8iFMAcfqfnP3JtVOgVLNlbyELVu7i7YwsCorLCQ8zeH0nb6nhz90Dg1l9u7NSjTHPA/dba8uNMR2Bp4AkQElEQkdRrrOaPP0lyP0WWrV11nWMuQ56jdEMKzlBUWkF763NZuGq3WTsyqdVeBgXDevOleMSyCko4b7F60/o0vL37oHBrL5JZDTwPLCqcnD9dzi1s5L8G5ZIAPh8sO1TJ3Fsfh985RA3DmY/6awmjwr+EhPSdKy1rPUUsHDVLpZkZlNU5mVAbFvunzGES8fE0bnN8c3BTJgJ6O6BwazeOxsaY2KAb4BhwLPW2lsCEVigaGfDFqjA4ywGzHgFCnY5NaxGXQljrtVmT3KSgiPlLM7MYuGq3Wzac4joyDAuGdmLK8fFMyahU9AWQgw0v+xsaIwZDbyCs7jwPuAvxpjXgDuqLEAUcd/RBYHpLzkLAq0P+k+GC38Pg2dAROjUJpLAs9aycvsBFq7azQfr9lBa4WNE7w78Yc5wZo3uRfvo4Ni7IxjVtztrGfBLa+0z4OyzDjwOrMNZMyLirtytkHF0QeB+aNcLJtwNSVdDp75uRydBJvdwKYvSPLy+ajfbcotoFxXBZSlxXDE2geG9O7gdXkiobxIZa63ddvSBtbYIuMkYM8u/YYnUQ3kxbHynckHgV86CwMTpziD5GVO0IFBO4PVZvtyay8KVu/ho4z4qfJaUPp24Y/IAZozoSUyr0CuC6Kb6rhPZZoy5ELgS6GatnWmMSQEOByQ6kdPZs7ZyQeAbzt7knfvDBb+DUVdBu+5uRydBZk9BMW+s8vDG6t1k5RfTqXUkN5zdlyvGxTMgVnvZN1R9x0TuBH4KPAPMqzxcjNOldbZ/QxOpQUkBrKtcELgn8/iCwOTrnQWBLXTQU2pW7vXxyeYcFq7cxeff7sdn4dwBXbn34sFcOLQ7URFqdTRWfdv5dwFTrLU7jDG/rDy2Gac8vEhgHF0QmPaisyCwohi6D4fpj8DIyyCmk9sRSpDZmVfEwlW7SU3zsL+wlO7to7hj0gAuHxtPfOfWbofXrNQ3ibQDdlf+/ejc4EigzG8RiRx10oLAds4+HWOug15JanXICY6WIVm4cjcrtuURZuD8wbFcMTaBSYndiAhX2ZpAqG8SWQ78CniwyrH/AD71W0TSsvl8sO2TygWBHzgLAuPHOwsCh82FVm3cjlCCzLf7jpchyT9STnznGO6+aBDzk+Pp0aF5lV0PRvVNIncC7xpjfgS0M8ZsAQqBS/wembQsNS0IHH8rJF0LsYPdjk6CTFFpBe+v3cOCVbvI2JVPZLjhomE9uHJsAmef0SWkdwoMNfWdnbXHGDMWGIuzNe5uYKW11heI4KSZ85bDln8eXxCIdRYEXvRfkHixFgTKCay1rMsqYMHK3by7JpvDpRWnLEMiTafeE+itUydlZeWXSP3lfuckjjULji8IPO8eLQiUGhUUl/NOZhYLVp5YhuSKsfEk92m5ZUiChVZhSdMoOwKbljgzrHZ9DWERlTsEXg8DpkCYplrKcdZaVu04yMKVu3i/sgzJ8N7t+e85w5mtMiRBRUlEAsdayEqHjJdh/VtaECi1yjtcyqJ0DwtX7WbbfpUhCQVKIuJ/h3OcbWUzXoH9myEixlkQmHQN9D1XU3NbuMUZWSeUTb/7wkF0aRfFwlVOGZJyr1OG5Pb5ZzBjZE9at9KvqWDW4H8dY8xya+15/gxGQpi3HL770Jlh9d1S8FU4e3XM/CsMuxSi27sdoQSBxRlZ3PvWumMbOGXlF/PzN9dggU6tI7n+rL5cPjaegd1VhiRUNCbFn+O3KCR05WxyWhxrX3cGydt2h7N+DKOvhm4qZCAn+tO/Np+wAyBwLIH8+74pKkMSgtROlPorKYD1i5zkkZV2fJA86VoYcIGq5soJrLWszzrEm2m7yS4oqfGc/CPlSiAhSv/bpW58Ptix3Omu2rQEKkogdihM/SOM+AG07eZ2hBJk9heW8k5mFm+u9rBlXyGtIsKIiQw/qSUC0KtjjAsRij8oicjpHdzpbPCU+Zqzkjy6gzNAPvpq1a+Sk5RVOFVzU9N28+mW/Xh9ltHxHXlw7nAuGdmLTzfnnDAmAhATGc49U9X1GaqURORk5cWw6V1nau725YCB/pPgggecrWUj9alRTrQhu4DUNA/vZGZzoKiM2HZR3DyhH5clx52wV8ecpN4AJ8zOumdq4rHjEnoak0T0EbQ5qWlNR8c+MPnXMOpK6Kjdj+VEeYdLeSczm9Q0Dxv3HKJVeBgXDu3O/JQ4JgzoesqquXOSeitpNCONSSKf+y2KRjDGzAFmAO2BZ621H7obUYg53ZqOPudAmMpny3HlXh+fbdlPatpuPtmcQ7nXMjKuA/81exizRvWiY2vVr2ppGpxErLWTG/vmxpjncCoA51hrh1c5Pg34KxAOPGOtffg0cSwGFhtjOgGPAkoitalxTcfYyjUdc51xD5Eqtuwt5M3Vu1mcmUXu4TK6tm3FDWf3ZV5yHIN7aA1QS+b2mMgLwBPAS0cPGGPCgSeBCwEPsMoYswQnoTxU7fU3WmtzKv9+f+Xr5FSqr+loEwtn3uG0OrSmQ6rJP1J2rLtqXVYBEWGGKUNiuSw5nomJ3YjUJk+Cy0nEWrvcGNO32uFxwFZr7TYAY8xCYLa19iFq2LfEOCU8Hwb+aa1ND3DIoeeUazquqVzToUJ2clyF18cX3+XyZtpuPt6YQ5nXx9Ce7fntJUOZPboXXdqqPL+cqDFlTwxwvbX2Bf+FA0Bvjm/BC05rZPxpzr8TuADoYIwZYK39e/UTjDG3ALcAJCQk+DHUIKU1HVJPW3MKeTPNw9vpWeQUltK5TSuuPjOB+clxDOul7k05tcaMiVhjzEycLinXWGsfBx6v5ZyngacBUlJS7OnODWnV13REdXDWcyRdozUdcpKC4nLeXZPNm2ke1uzOJzzMMDkxlvnJcZw/OJZWEequkto1tjsr0hjzGbAa8AFYa/+zkdfMAqrOJ42rPCY1KTvirOnIfEVrOqRWXp/ly625pKZ5WLphL2UVPhK7t+P+GUOYPbo33dqpu0rqp7FJ5M/VHvvjU/4qYKAxph9O8rgCuMoP120+rIXdK53Esf5tKCt01nRMug9GXwkdW0CXnZxS9VLr90xNZGRcB1LTPLyVnsXeQyV0iInkirHxXJYcz/De7bU7oDRYg5KIMeZWa+1TOAPd1RPH8npcZwEwCehqjPEAD1hrnzXG/ARYijMj6zlr7YaGxNnsFGTB2oVOd1XeVohsA8PmwOirIOFsremQGkut/+yNTKyFMAMTB3XjtzOHMmVIrAoeil80tCXy78o/32vMm1trrzzF8Q+ADxpz7WajvBg2v+8kjm2fgvU5iwDP/bmzKDCqrdsRShD509IaSq1baB8dwcc/n0hs+2iXIpPmqkFJxFq7pvKv64EDlYPsBujst8haMmud6biZr8K6RU4Jkg7xMOFup7uqc3+3I5Qgs23/YRale8jOr7nUemFJhRKIBERjx0TesNZOgWOztd4ApjQ+rBaqcC+sqeyuyt1SWYJkljPDqu8EdVfJCQqKy3l/7R5S03aTviufMANREWGUVvhOOlel1iVQGj07q9pjFc6pr4pS2PJPp9Wx9WOnuyp+PMx83BnvUAkSqeLo7KpFlbOrSit8DIxty73TBzMnqTcrvs9TqXVpUo1NImuNMX/FKcY4EVjb+JBaAGthT6bT4lj3JhQfhHa94Jy7nFZH1wFuRyhBZmuO0131dpXZVT9IiWd+chwj4zocm12lUuvS1Iy1tc/KNcbcA6QDGdbaA9WemwkMATZZa98NSJR+lJKSYlevXu3Omx/eD+vecFaS52yA8CgYcokzu6r/ZAjTbBk5ruBIOe+udWpXZVYuBpw4qBvzk+M0u0qanDEmzVqbUv14XVsiU4FfAp0qp+Km46zneKcycQR98nCNtxy+Xeq0Oo5WzO2dDDMeg+GXQkwntyOUIFLh9fFF5WLAjzbuO7YY8NcXD2F2Ui9i22lwXIJLnZKItfYCAGNMHyAJGAOcB/ymssLuD621RwIWZSjau94Z51j7BhzJhbbdnYq5o6+G2MFuRydB5rt9haSmH69d1bF1JFeNO1q7SosBJXjVa0zEWrsT2AksBjDGdAVewynDfp+/gws5RXmwPtWpmLt3LYS3gsTpTuI4YwqEu115X4JJ/pEy3l3jdFet8RRUqV3Vm8mD1V0loaFRv9WstbmVq8s/oKUmEW+FM6sq81VnlpWvHHqOhumPwIj50FpLZ+S4Cq+P5d/tJzXNc6zU+uAeql0locsfH413AT39cJ3QkrPZqV215nUoyoHWXWHcLc4geY/htb9eWpQtewtJTdvN2xnZ5B5WqXVpPuqURIwx+UBG5Vd65Z+brLU+4Grg+0AFGFSKD1Zu8PQqZKcf3+Bp9FUw8CJt8CQnOFBUxpLMLBalZx3bGfD8wU6p9UmJKrUuzUNdWyIzcAbUk4CfAcMAa4wpBqKAHwQmvCDy8e9hxZPgLYXuw2HqQzDiMm3wJCco9/r4bMt+FqV5WLZ5H+Vey7Be7Xlg5lBmjdLOgNL81HV21lfAV0cfG2MigKFALLDOWrsvMOEFkY4JkHwDJF0NPUZqg6cWrKZS64k92pGa5mFxRhZ5RWV0bduK687qy7wxcQzt1d7tkEUCpk6LDZsTVxcbSsirXmodwODshxAZbpgyuDvzk+OYmNiNyHB1V0nz0eDFhsaY6+r5XpnWWpU/kWbpT/+qodQ60CEmks/unkSnNiofJy1LXbqzJlO3HQuPfiDLRzW0pBmx1rIh+xCpaR6yC2outX6ouFwJRFqkuiSRHTjJoS6DAEeTiEjIyyks4Z2MbBale9i8t5BW4WHERIZRXK5S6yJH1XV2Vl1HkTXaLCGttMLLsk05pKZ5+Pzb/Xh9ltHxHfnvOcOZObInn23Zr1LrIlXUmkSstb9vikBE3GKtZY2ngEVpHpasyaaguJzu7aO45bz+zBsTx4DY41sQq9S6yIlUzElarL0FJbydkcWidA9bcw4TFRHG1GE9mJ8cxzkDuhIeVnPDek5SbyUNkUpKItKilJR7+XDjPlLTPHz53X58FlL6dOKhS0cwY2RP2ker6oBIfSiJSLNnrSV910FS07J4b202hSUV9OoQzY8nD+DSMXH069rG7RBFQpaSiDRbWfnFvJ3uYVF6Fttzi4iJDGf6cKe76sz+XQg7RXeViNSdkog0K0fKKli6YS+paR6+/j4Pa2F8v87cPukMLh7Rk7ZRuuVF/En/oyTkWWtZuf0Ai9I9vL92D0VlXuI7x/DTKQO5NCmOhC6t3Q5RpNlSEpGQtfvAEd5Kd2ZX7TpwhDatwrl4RE/mJ8cxtm9ndVeJNAElEQkpRaUVfLBuD6lpHr7ZfgBj4OwzunDXBQOZNrwHrVvplhZpSvofJ0HP57P8e1seqeke/rluL8XlXvp2ac3dFw1i7pg4eqvkiIhrlEQkaO3ILWJRuoe30rPIyi+mXVQEc5J6MT85jjEJnTDa00XEdUoiElQOlZTzwVqnu2r1zoMYAxMGduM/pyUydVgPoiPD3Q5RRKpQEhHXeX2Wr7bmsijdw7/W76W0wscZ3drwy2mDmZvUmx4dot0OUUROQUlEAq6m7WTnJPVma85hFqV7eDs9i72HSmgfHcFlKXHMT45nVFwHdVeJhABtjysBVdN2spHhhp4dotl1oJjwMMPEQd2YNyaOKUNi1V0lEqQavD2uSGM8snTLSdvJlnst2fkl/PriIcxO6kVsO3VXiYQqJREJmC17C8nKL67xOa/P8qPz+jdxRCLib0oi4lcHi8pYsiab1DQP67IKTnmetpMVaR6URKTRyr0+Ptuyn0VpHpZt3ke51zK0Z3t+e8lQIiMMf3x/s7aTFWmmlESkwTZmHyI1zcM7mVnkFZXRtW0rrjurL/PGxDG0V/tj57WLitR2siLNlJKI1Evu4VLeyXS6qzbtOURkuOGCId2ZNyaOiYndiAwPO+k12k5WpPlSEpFalVX4+GTzPlLTsvhsSw4VPsvIuA781+xhzBzZi05tWrkdooi4JOSTiDGmDfA58Dtr7Xtux9NcWGtZn3WI1LTdvLMmm/wj5cS2i+Kmc/sxLzmOQd3buR2iiAQB15KIMeY54BIgx1o7vMrxacBfgXDgGWvtw7Vc6pfAGwELtIXJOVTC4swsUtM8fLvvMK0iwrhoaHfmJccxYUBXImrorhKRlsvNlsgLwBPAS0cPGGPCgSeBCwEPsMoYswQnoTxU7fU3AqOAjYBWqzVCSbmXjzftIzXNw/Jv9+OzMCahIw/OHc4lI3rRoXWk2yGKSJByLYlYa5cbY/pWOzwO2Gqt3QZgjFkIzLbWPoTTajmBMWYS0AYYChQbYz6w1vpqOO8W4BaAhIQEP34XoctaS8bufBaleXh3TTaHSiro2SGa2yedwaVj4jijW1u3QxSREBBsYyK9gd1VHnuA8ac62Vr7awBjzA1Abk0JpPK8p4Gnwamd5a9gQ9GeguJjW8pu219EdGQY04b1YH5yPGed0YVwbSkrIvUQbEmkQay1L7gdQzArLvPy4ca9pKZ5+HJrLtbCuL6dufW8/lw8oiftotVdJSINE2xJJAuIr/I4rvKYnEZNpdZnj+7F6p0HWZTm4b21ezhcWkHvjjHcef5A5o3pTZ8ubdwOW0SagWBLIquAgcaYfjjJ4wrgKndDCm7VS61n5Rdz95tr+MN7G8gtKqd1q3CmD+/J/OQ4xvfrTJi6q0TEj9yc4rsAmAR0NcZ4gAestc8aY34CLMWZkfWctXaDWzGGgppKrVf4LIdKvTx62SimD+9Bm6hg+6wgIs2Fm7OzrjzF8Q+AD5o4nJDj81m+2X7glKXWyyt8zE+Oa+KoRKSl0UfUELMzr4hF6Vm8le7Bc7AYA9Q03Uyl1kWkKSiJhIDCknI+WLeHRWlZrNxxAGPg3AFduWdqIqXlPh5YskGl1kXEFUoiQcrrs6z4Po/UtN38a8NeSsp99O/WhnumJjI3qfcJLY1WEWEqtS4irlASCTLb9h9mUbqHt9Kz2FNQQrvoCOaNiWNechxJ8R0x5uTZVSq1LiJuURIJAgXF5by31tmjI2NXPmEGzhvUjV/PGMIFQ7oTHRnudogiIjVSEnFJhdfHF1tzWZTm4cON+yir8DGoe1vunT6YuUm9iW2vmpIiEvyURJrYt/sKWZTm4e2MLHIKS+nYOpIrx8YzLzmOEb071NhdJSISrJREmsDBojLereyuWuspIDzMMDkxlvnJvZk8OJaoCHVXiUhoUhIJkHKvj8+37GdRuoePN+2j3GsZ0rM9v7lkKLNH96Jr2yi3QxQRaTQlET/bmH2IReke3snMIvdwGV3atOLaM/syL7k3w3p1cDs8ERG/UhLxg9zDpbyTmc2iNA8b9xwiMtwwZXB35ifHMTGxG5HaUlZEmiklkTqoqdT6xSN68snmHFLTPHy2JYcKn2VkXAd+P2sYs0b1olObVm6HLSIScMbalrXRX0pKil29enWdz69eah0gPMwQFW44Uu6jW7soLk3qzbzkOAZ1bxeIkEVEXGeMSbPWplQ/rpZILWoqte71WWy44fkfjmXCgK5EqLtKRFooJZFaZJ+i1HpJuY/JibFNHI2ISHDRR+hanKqkukqti4goidTqnqmJxFSrXaVS6yIiDnVn1eJodVyVWhcROZmSSB2o1LqISM3UnSUiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg3W4mpnGWP2AzuBDkBBPV9en9fUdm5jnj/Vc12B3DpF546G/Myb8tr1vYY/74faztH90PTXDuTviMbeD6d7PlD3Qx9rbbeTjlprW+QX8HQgX1PbuY15/lTPAavd/rn6+2felNeu7zX8eT809N9c90Pw3A/1eU1j74da/t2b9H5oyd1Z7wb4NbWd25jnGxJ7MAhk3P64dn2v4c/7obZzdD80/bUD+TuisffD6Z5v0vuhxXVnNWfGmNW2hlLN0jLpfpCqAnU/tOSWSHP0tNsBSFDR/SBVBeR+UEtEREQaTC0RERFpMCURERFpMCURERFpMCWRFsIYM8cY8w9jzOvGmIvcjkfcZYzpb4x51hiT6nYs4g5jTBtjzIuVvxeubuh1lERCgDHmOWNMjjFmfbXj04wxW4wxW40xvzrdNay1i621PwJuAy4PZLwSWH66H7ZZa28KbKTS1Op5b1wKpFb+XpjV0PdUEgkNLwDTqh4wxoQDTwLTgaHAlcaYocaYEcaY96p9xVZ56f2Vr5PQ9QL+ux+keXmBOt4bQBywu/I0b0PfUDsbhgBr7XJjTN9qh8cBW6212wCMMQuB2dbah4BLql/DGGOAh4F/WmvTAxyyBJA/7gdpnupzbwAenESSSSMaFGqJhK7eHP8UAc4Ncbo9fO8ELgDmG2NuC2Rg4op63Q/GmC7GmL8DScaYewMdnLjqVPfGW8A8Y8zfaESpFLVEWghr7ePA427HIcHBWpuHMz4mLZS1tgj4YWOvo5ZI6MoC4qs8jqs8Ji2T7gc5lYDeG0oioWsVMNAY088Y0wq4AljickziHt0PcioBvTeUREKAMWYBsAJINMZ4jDE3WWsrgJ8AS4FNwBvW2g1uxilNQ/eDnIob94YKMIqISIOpJSIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCIiIg2mJCLiMmPMJmPMYWNMWeXX4cqvIW7HJlIblT0RCRLGmGeBbdbaB92ORaSu1BIRCR4jgfW1niUSRJRERIKAMSYMZ/9rJREJKUoiIsEhAef/4za3AxGpDyURkeDQHigCWrkdiEh9KImIBIdNwBrgoDFmsNvBiNSVZmeJiEiDqSUiIiINpiQiIiINpiQiIiINpiQiIiINpiQiIiINpiQiIiINpiQiIiINpiQiIiINpiQiIiIN9v8BjfXWng8Pia0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "taus = np.array([0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1])\n",
    "errors = []\n",
    "\n",
    "for tau in taus:\n",
    "    errors.append(calculate_total_error(tau, order=2)) # 通过 order=2 指定 Suzuki 分解的阶数    \n",
    "\n",
    "plt.loglog(taus, errors, 'o-', label='error')\n",
    "plt.loglog(taus, (2 * taus * 3 )**3 / 3 * (1/taus) * np.exp(3 * taus), '-', label='bound') # 按照 (12) 计算的二阶误差上界\n",
    "plt.legend()\n",
    "plt.xlabel(r'$\\tau$', fontsize=12)\n",
    "plt.ylabel(r'$\\Vert U_{\\rm cir} - \\exp(-iHt) \\Vert$', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c468fd4",
   "metadata": {},
   "source": [
    "可以看到，实际计算出的模拟误差都小于其理论上界，说明这样的结果是符合预期的。\n",
    "\n",
    "## 小结\n",
    "\n",
    "本教程介绍了 Paddle Quantum 中构建模拟时间演化电路的功能，并对其背后的理论基础——Suzuki product formula 进行了介绍。利用该功能，用户可以通过搭建对于自定义哈密顿量的任意阶 product formula 电路来模拟不同物理系统的时间演化过程。\n",
    "\n",
    "量子模拟本身是一个比较宽泛的话题，其应用也十分广泛：凝聚态物理中的多体局域化、时间晶体、高温超导、拓扑序的研究；量子化学中的分子动力学模拟、反应模拟；高能物理中的场论模拟；乃至核物理和宇宙学中的相关研究。本教程中介绍的 Suzuki product formula 与基于通用量子计算机的数字量子模拟只是量子模拟的一部分，许多不基于通用量子计算机的量子模拟器，例如在冷原子、超导、离子阱以及光子等平台上做的模拟量子模拟（analogue quantum simulation）也是一个十分重要的方向。鉴于篇幅，本教程中无法进一步对这些内容逐一展开介绍，感兴趣的读者可以阅读 14 年的这篇综述文章 [8]。对于一些更新的结果，也可以参考这篇中文综述 [9]。\n",
    "\n",
    "在后续的教程 [模拟一维海森堡链的自旋动力学](./SimulateHeisenberg_CN.ipynb) 中，我们以凝聚态物理中的自旋模型为例，进一步地展示了如何利用本文中介绍的内容来进行量子多体模型的动力学模拟。同时，该教程也将介绍如何搭建不基于 Suzuki 分解的自定义时间演化电路。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01705a8d",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 参考资料\n",
    "\n",
    "[1] Feynman, R. P. \"Simulating physics with computers.\" International Journal of Theoretical Physics 21.6 (1982).\n",
    " \n",
    "[2] Lloyd, Seth. \"Universal quantum simulators.\" [Science (1996): 1073-1078](https://www.jstor.org/stable/2899535).\n",
    "\n",
    "[3] Childs, Andrew M., et al. \"Toward the first quantum simulation with quantum speedup.\" [Proceedings of the National Academy of Sciences 115.38 (2018): 9456-9461](https://www.pnas.org/content/115/38/9456.short).\n",
    "\n",
    "[4] Nielsen, Michael A., and Isaac Chuang. \"Quantum computation and quantum information.\" (2002): 558-559.\n",
    "\n",
    "[5] Tran, Minh C., et al. \"Destructive error interference in product-formula lattice simulation.\" [Physical Review Letters 124.22 (2020): 220502](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.220502).\n",
    "\n",
    "[6] Childs, Andrew M., and Yuan Su. \"Nearly optimal lattice simulation by product formulas.\" [Physical Review Letters 123.5 (2019): 050503](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.050503).\n",
    "\n",
    "[7] Campbell, Earl. \"Random compiler for fast Hamiltonian simulation.\" [Physical Review Letters 123.7 (2019): 070503](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.070503).\n",
    "\n",
    "[8] Georgescu, Iulia M., Sahel Ashhab, and Franco Nori. \"Quantum simulation.\" [Reviews of Modern Physics 86.1 (2014): 153](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.153).\n",
    "\n",
    "[9] 范桁. \"量子计算与量子模拟.\" [物理学报 67.12(2018):16-25](http://wulixb.iphy.ac.cn/article/id/72211)."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3b61f83e8397e1c9fcea57a3d9915794102e67724879b24295f8014f41a14d85"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
